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Journal of the 
ENGINEERING MECHANICS DIVISION 


Proceedings of the American Society of Civil Engineers 


MATRIX ANALYSIS OF BEAMS 


Ray W. Clough,! A.M. ASCE 
(Proc. Paper 1494) 


SUMMARY 


Matrix formulation of the general theory of structures has been the sub- 
ject of several recent papers. In the present paper, the application of a 
matrix procedure to the analysis of structures consisting of flexural mem- 
bers is presented. The procedure may be applied to non-uniform as well as 
uniform beams, and permits consideration of shear distortion effects. Four 
examples demonstrate the application of the method to various types of 
structures. 


Notation 


Symbols are defined where they are first introduced in the text, and the 
most important ones are also listed here for convenience. Underlining indi- 
cates that the symbol represents a matrix. A “prime” indicates a transposed 


matrix. 

fn = flexibility of element “n” 

f = flexibility of complete system of elements 

bo = element forces in the primary system due to unit 
external loads 

by = element forces in the primary system due to unit 
redundants 

D = deflection of redundants due to unit redundants 


Note: Discussion open until June 1, 1958. A postponement of this closing date can be 
obtained by writing to the ASCE Manager of Technical Publications. Paper 1494 is 
part of the copyrighted Journal of the Engineering Mechanics Division, Proceedings 
of the American Society of Civil Engineers, Vol. 84, No. EM1, January, 1958. 


1. Associate Prof. of Civ. Eng., Univ. of California, Berkeley, Calif. 
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values of redundants due to unit external loads 


flexibility matrix of the complete structure 
stress matrix of the complete structure 
external load system 

length of beam segment 

rotation of end of beam segment 

linearized rotation of end of beam segment 
moment at end of beam segment 

linearized moment 

beam flexural stiffness 


modules of rigidity 


effective shear area of beam 


bending flexibility factor 


= shear flexibility factor 


elastic support flexibility factor 


spring constant of elastic support 
force acting in elastic support 
external force applied to structure 
external moment applied to structure 


linearized external moment 


superscripts referring to left end, center, and right 
end of a beam segment respectively 


subscript indicating the number of the element 
considered 


= subscript indicating reference quantity 


(It will be noted that the subscripts in the symbols by and by are excep- 
tions to the above usage, but this should not cause any confusion.) 


INTRODUCTION 


At least two general approaches are available for the systematic analysis 
of beam flexibility and stresses. One of these is essentially an extension of 
the Holzer method of evaluating the critical frequencies of shaft vibrations. 
The initial developments in this method were made by Myklestad,(1) while 
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Thomson(2) and Marguerre(3) have concerned themselves with matrix formu- 
lations of procedure. This approach is well suited to the study of vibrations 
of beam systems, but is somewhat awkward to apply to static analyses. 

On the other hand, the standard procedures of structural analysis, such as 
the virtual work (or dummy unit load) method provide very convenient means 
for evaluating stresses or deflections under static conditions. The method 
presented herein may be considered as a matrix systemizing of the virtual 
work procedure. The numerical operations which are involved correspond 
exactly to those which would be performed in the standard virtual work solu- 
tion. The advantage of the matrix formulation lies in the fact that all prob- 
lems may be treated by an efficient standardized procedure. Moreover, the 
mathematical operations necessary after the original matrices have been 
established are of such routine nature that they can be carried out by persons 
having no knowledge of structural analysis, or preferably by an electronic 
digital computer. 

The direct results of this matrix procedure are two influence coefficient 
matrices: the first of these, the stress matrix, represents the forces in each 
of the members of the structure due to successive applications of unit values 
of the external loads; while the second, the flexibility matrix, represents the 
deflections at the points of loading due to unit values of the external loads. 

. When these two matrices have been evaluated, it is a simple operation to ob- 
tain stresses or deflections due to any system of externally applied concen- 
trated loads. (As is true of influence coefficient procedures in general, this 
method does not, without modification, make possible the analysis for distri- 
buted loads. Usually the simplest manner of handling a distributed load is to 
replace it by an equivalent system of concentrated loads applied at points for 
which influence coefficients have been obtained.) 

The matrix formulation of structural theory which forms the basis for 
this paper was pioneered by Langefors(4) and recently was given a very 
thorough and effective treatment by Argyris.(5) The purpose of this paper is 
to demonstrate the application of the Argyris matrix equations to the analysis 
of stresses and deflections of various types of beam systems. A convenient 
approximate procedure for taking account of non-uniform beam stiffness is 
also presented and shown to give good results. 


General Outline of the Method 


In this presentation, the Argyris notation, as summarized by Hunt(6) will 
be used to represent the various matrices involved in the analysis. The 
principal matrix equations will be restated here for purposes of clarity. 

The procedure involves first establishing a flexibility matrix f and two 
force transformation matrices bg and bj. The flexibility matrix is diagonal 


matrix of submatrices f1, fg, f3, . . . fp, thus 
fo 
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where the submatrices fj, etc., represent the flexibility of the “n” elements 
into which the structure is considered to be divided for purposes of analysis. 

To establish the force transformation matrices, it is first necessary to 
make the system statically determinate in the usual way, by selecting a set 
of redundant forces (one for each degree of static in determinacy). The force 
transformation matrix b, then represents the static equilibrium relationships 
between the external loads applied to the structure and the forces developed 
in each element, the redundants being assigned arbitrary values (usually zero). 
Similarly, b} represents the equilibrium relationship between the redundant 
forces in the structure and the element forces. The by matrix exists only 
for statically indeterminate structures, of course. Generally several differ- 
ent by matrices could be set up for a given structure depending upon the 
particular set of forces which are selected as redundants. As is well known, 
the conditioning of the matrix which must be inverted in the solution of the 
structure will be influenced greatly by the choice of the redundants, but such 
considerations are not a part of the present discussion. 

When the flexibility of the elements and force transformation matrices 


have been established, the analysis of the structure is carried out by a 
series of matrix operations as follows: 


i) Evaluate the displacements corresponding to the redundants* 


D fb 
2) Compute the values of the redundants due to application of unit external " 


loads 
(A-2) 


3) Compute the stress matrix, i.e. the total forces in all the elements of 
the structure due to unit external loads 


b= & * 


4) Compute the flexibility matrix of the structure, i.e. the displacements 


corresponding to the external loads due to application of unit external 
loads 


(A-4) 


5) Check the structural continuity at the redundants when unit external 
loads are applied 


= (A-5) 
After the stress and flexibility matrices have been obtained from Eqs. A-3 


and A-4, the stresses S and deflections r due to any system of applied loads 
R may be calculated from following equations 


S- bR (A-6) 


*The structure of this equation is discussed in the Appendix. 
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It will be noted that the procedure involves only matrix multiplications and 
additions except for the inversion of the matrix D. This operation is 
equivalent! to the solution of the set of simultaneous equations obtained from 
the conditions of continuity imposed at the redundants. There will be as 
many conditions as there are redundants, hence the analysis of a structure 
indeterminate to the degree “m” will require the inversion of a matrix of 
order “m.” Generally this inversion will cause little difficulty if the calcula- 
fioh is carried out on an automatic digital computer, but if the analysis is to 
be done by desk calculator the inversion may represent a large part of the 
total work to be done. In this case a reduction in the size of the matrix to be 
inverted may be very desirable. Langefors has pointed out(7) that such a 
reduction may be accomplished by analyzing parts of the structure first, and 
then combining these parts in a subsequent operation. (G. Kron has also 
made use of this concept in his “tearing” procedure, described in References 
10 and 11.) In effect, the flexibility matrix of each subsection is obtained by 
following the procedure outlined in Eqs. A-1 to A-5. Then the subsections 
are joined by repeating the procedure, the subsections now being considered 
as structural elements. In order to carry out this joining operation, the 
forces acting at the points of contact of the subsections must be included in 
their flexibility matrices. In other words, the interconnection forces between 
the subsections must be considered as part of the external force systems 
acting on the subsections, even though these forces may not be required in 
the combined flexibility matrix. An analysis of this type is demonstrated in 
Example IV of this paper. 


Flexibility of the Beam Segments 


Since this presentation is limited to the analysis of beam type structures, 
only the flexibility of beam segments (and possible elastic supports) need be 
considered here. For the purposes of this discussion, it is convenient to 
define the flexibility of a beam segment in terms of the rotations at the ends 
of the segment developed by unit moments applied at the ends (see Fig. 1). 
Rotations and moments are considered positive in the directions shown in 
the figure. In matrix notation, the rotations of a beam segment in terms of 
the moments, are as shown in Eqs. 1, in which y is the flexibility matrix 
of this beam segment. When flexural elements are combined with linear 
elements, such as elastic supports, it is somewhat more convenient to ex- 
press the flexibility in terms of linear forces and displacements; accordingly 


Fi‘ 


t, (la) 


1. Further discussion of this point is given under “Discussion” at the end 
of the paper. 
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Fig. 1 Typical Beam Segment “rn 


Eq. 1 will be written as shown in Eqs. 2. In Eqs. 2b, 1, is any convenient 


reference length. 
oh git ‘ 


(2a) 
where % 2 


™.= Ma/ee 


or 


(2b) 


For calculation purposes, the flexibility f, may be expressed as the sum 
of two matrices: flexibility due to bending {Bp and flexibility due to shear 
fsn , thus 

The elements of these matrices represent the rotation at the ends of the seg- 
ment due to unit moments applied at the ends of the segment. The bending 
flexibility terms may be calculated from the expressions of Eqs. 4. The 
values of these terms will depend, of course, upon the manner of variation 
of the stiffness El(x). Various assumptions may be made regarding how 
this quantity varies, with corresponding differences in the results; however, 
only two cases will be considered here. If an average value for the stiffness 
ET, is assumed constant through the length of the segment, the flexibility 
matrix will take the form of Eqs. 5. The symbols l, and Elo in Eq. 5a 
represent any convenient reference values of length and stiffness. 
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On the other hand, if it is assumed that the flexibility a varies para- 


bolically through the length of the segment the flexibility matrix is given by 
Eq. 6. EI,', EI,© and EI, being 


EL, |(O.lak + +O.l0k) (-O.lan + L208 


where @, EI. 


ET, 
Ele 
Els 
Ele 


@. 


values of the stiffness of the beam segment at the left end, middle, and right 
end respectively, and lo and Elo being reference values as before. Eq. 6 
gives much better results than Eq. 5 if there is an appreciable variation of 
stiffness within the beam segment, and for such cases the extra effort in- 
volved in its evaluation is generally worthwhile. 

Elements of the shear flexibility matrix may be evaluated similarly as 
shown in Eq. 7 in which G is the modulus of rigidity of the material and A 
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represents the effective shear area of the cross-section, (see for example, 
W.M. Fife and J. B. Wilbur, “Theory of Statistically Indeterminate Structures,” 
McGraw-Hill, New York, 1937, p. 29). For a rectangular beam section Al is 
5/6 of the total area, while for flanged beams it may be taken equal to the web 
area. As is apparent from Eq. 7, the shear sens weed matrix depends only 


upon the average value of the shear flexibility, GAL , and may be expressed 
as shown in Eq. 8. 


-l 
f. = 


= GEI. 
SA 


(Ba) 


where 


For a parabolic variation of the shear flexibility, the average value within 


the segment is given by Eq. 9 in which the superscripts have the same mean- 
ing as before. 


cA, | (a) * (am) + | (9) 


Eqs. 5 and 8 may be combined to give the total flexibility of the beam seg- 
ment. Where bending flexibility may be represented by an average value, the 
combined flexibility matrix takes the form shown in Eq. 10a while a parabolic 
variation of the bending flexibility yields the matrix of Eq. 10b. 


C Xn Gn) ( 


+ 0.8x5 - Bn) + 1.208 + + Bn) 


Eq. 10a or 10b provides the only type of element flexibility matrix needed 
for the analysis of structures composed of beam or girder elements. In 
cases where the beams are mounted on elastic supports, however, the flexi- 
bility of the supports must also be considered. A linear elastic support of 
any type may be represented by a spring as shown in Fig. 2. The deflections 
of such an element may be represented as shown in the following equation, 


Sn = fy, Pr where (b) 


k, being the spring constant. Since only one type of displacement is con- 
templated, the flexibility matrix for such an element consists of a single 
term. In order to facilate the combination of elastic supports into a beam 


we 
is 
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n 


Spring Constant 


n 


Fig. 2. Elastic Support “n* 


system, it is convenient to express the spring flexibility in the form given in 
Eq. 11, Ely 


3 
Ey, Tn (11) 
— GEL 
where GEI./e2 
= (11a) 


and |, being reference values as before. The analysis of any type of beam 
system may now be carried out using the element flexibilities represented by 
Eq. 10a or 10b and Eq. 11. 


Examples of Beam Analysis 


I. Uniform Simple Beam 


As an introductory example the deflections and rotations will be evaluated 
at a series of equally spaced points across the span of a uniform beam, due to 
vertical forces and moments applied at these points, as shown in Fig. 3. 


EI 


Fig. 3 Uniform Simple Beam 


ike 
| 
q 
x 
al 
a P p Pp Pp 
2 3 + 5 
q 
l 
. 
‘tol 
if 
late 


1494-10 EM1 January, 1958 


Flexibility Matrix: 
In this uniform structure, all beam segments have the same flexibility, as 
shown in Eq. 12a (shear deformations being neglected). 


The complete flexibiljty matrix for all beams then becomes 


e? 001 2 (12b) 
« 


Force Transformation: 
For this statically determinate system, the redundant force matrix bj 
does not exist. Thus, the matrix b, represents the complete relationship be- 


tween the applied loads and the member forces in terms of which the flexibili- 
ty matrix has been expressed, as shown in Eq. 13a. 


™; | 000006000 00 0|fP, 
mi 5 4 39 214 -1| | 
5 43 21 5 § =1 91 
48 6 2 4 bh -2 P, 
1, 6 9 6 3 3 3 3 -3 | 
me 3 69 63 3 3 3 3 -3 -3| | A! (13a) 

24 6 86 2 2 2 2 | 
mi 12345 214 2 4 
mi: 000000 0 0 0 0 0 

where ™m, = Mn/y 

Jon = | 


In matrix symbols this equation may be written simply 


m= R 
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Solution: 
The flexibility matrix for a statically determinate system is given by the 
matrix operation 


be Ff (c) 


which is equivalent to Eq. A-4 since b, = b inthis case. Carrying out this 
operation yields the matrix shown in Eq. 14. Rows and Columns of this 
mairix are arranged in the sequence specified by the column matrix 

of external loads R in Eq. 13a. It may be used to evaluate the bending dis- 


placements and rotations induced by any system of applied loads and 
moments R , as shown by Eq. A-5. 


300 456 468 372 2041330 240 78 -48 -138 -192 -210 
! 
768 828 672 372/480 408 192 -60 -240 -348 -384 
972 828 4681486 432 270 0 -270 ~432 -486 
768 456,384 348 240 60 -192 -408 -480 
3001210 192 138 48 -78 -240 
j432 234 72 -54 -144 -198 -216 


(Symmetrical) 


Shear Flexibility: 

To demonstrate the effect of shear distortions on the beam flexibility, a 
second calculation was made in which the shear terms were incorporated into 
the flexibility matrix of the beam segments. For this example, it was as- 
sumed that the beam under consideration was a 24'' I 100# rolled section on 
a total span of 18 feet (1=3 ft.) for which the shear flexibility factor is given 
by 

GETZ = |.528 
GA’ 


making the total flexibility matrix of a single beam segment as shown in Eq. 
15. Since the geometry of the system has not been changed, the force trans- 


@ [3528 -osze 


formation matrix remains as shown in Eq. 13a. The flexibility of the beam 
was Calculated as before, with partial results as shown in Eq. 14a below 
(only the vertical displacement due to vertical loads being shown, to save 
space). 


q 
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zoe 90 =30 =120 | (24) 
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49267 49565 390.3 213-2 
841.3 883.0 708.7 390.3 


3 
L 1054.5 883-0 495-5 (ida) 
(26>*EI 
841-3 492-7 
345-8 


Comparing this with the corresponding portion of Eq. 14 it is evident that 
shear distortions have a significant effect on the deflections in this beam. 


Il. Tapered Cantilever Beam 


Application of the suggested approximate procedure for treating non- 
uniform beams will be demonstrated by analysis of the flexibility of the 
tapered cantilever beam shown in Fig. 4, for which an exact solution has 


been published by Bisplinghoff.(8) Only vertical loads will be considered and 
shear distortion will be neglected. 


Fig. 4 Tapered Cantilever 


Flexibility Matrix: 


Since the depth of the beam varies continously it was necessary to evalu- 
ate the relative stiffness at the ends and center of each segment. The stiff- 8 
ness at the support was taken as the reference value El,. The flexibility a 
matrix of each segment was then calculated according to Eq. 6. The flexibili- 
ty matrix of the complete set of beam segments is given in Eq. 16a. 
Force Transformation: 
For this statically determinate system the relationship between the applied 
loads and the moments in the segments is found easily to be as shown in 
Eq. 17a. At this point it may be noted that the size of the matrices (Eqs. 16a 
and 17a) may be reduced for this type of loading, due to the fact that the 
moments at adjacent ends of the segments are equal. (Mpj = My/,1). For 


this case the force transformation matrix of Eq. 17a can be reduced to the 
form in Eq. 17b. 


EI, | 
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2.27 1.30 0 0 0 0 0 0 
1.30 2.96 0 0 Oo 0 0 0 0 0 
0 0 3694 2631 0 0 0 0 0 
0 © 2631 0 0 0 0 0 0 
| o o 0 07670 4-73 Oo 0 
0 0 4073 11.64 0 0 0 o (| (16a) 
0 0 18.37 12.25 0 0 
0 0 12.25 32.68 0 
0 0 0 0 0 0 0 0 61.18 51.00 
0 0 0 0 51.60 164.76 
™: 1234 5 
0123 4 
™ 
0012 3| 
mi. |001 23 Po 
00012 
m: 0000 1 
00001 
or 
TH = be R 
4 
> 2 fe 6 4 2 3 P, (17b) 
Ms; 00001 


The flexibility matrix, Eq. 16a, may be put in the corresponding form by 
adding the rows (and columns) which correspond to the duplicated rows in 
Eq. 17a, i.e. 2 +3, 4+5, 647, and8 +9. The last row (and column) may 

be dropped because no moments are developed corresponding to these terms. 
The reduced flexibility matrix thus takes the form of Eq. 16b. When the 
calculations are to be done by automatic computers, there is generally no 
great advantage to be gained by this type of reduction, and the matrices may 
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2627 1-30 0 


1.30 6.90 231 0 0 
© 31 43:12 4.93 © (16b) 
0 0 4673 30.01 12.25 
0 oO 0 12.25 93.86 


as well be left in their original form. However, if the operations are to be 
performed with desk calculators, the reduced form is considerably more 

convenient to use. The remaining examples in this paper are expressed in 
the reduced form to save space. 


Results: 

Performing the required matrix multiplications (Eq. c) the flexibility 
matrix for this tapered beam was found to be as shown in Eq. 18a. For 
comparison purposes, the exact values given by Bisplinghoff, and values 


227 


5°84 9041 12.98 16.55 
21.18 38.83 56.48 74.13 


5 
F= Ea 85.99 137.88 189.77 (18a) 
(Symmetrical) 258.75 391.87 


obtained by assuming a constant average value of a for each section are 


given in Eqs. 18b and 18c below. (The average value for each segment was 
computed by a formula analogous to Eq. 9). 


From Bisplinghoff (8): 


2029 5686 9243 13.01 16.58 


2 21.21 38.85 56.48 74.12 
Fe = 86.02 137.93 189.90 (18b) 
P (Symmetrical) 259.10 392.40 


716.95 


Assuming Constant Segment Stiffness: 


2.6 


665 1004 1803 
2269 4167 60.5 7941 
91-9 146-8 201.8 


18c) 
GEI, (Symmetrical) 277.6 421.0 


| 

799.9 

= 
x 
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It is clear that the parabolic flexibility approximation is quite accurate in this 
case, whereas the assumption of uniform stiffness in each segment leads to 
Sizable errors. 

Ill. Beam Continous over Several Supports 


The application of this matrix method to a statically indeterminate system 
will be demonstrated first by the analysis of a uniform beam continous over 
four supports, as shown in Fig. 5. 


3 5 7 
2é 
Fig. 5 Uniform Continvous Beam 


Flexibility Matrix: 

Since the beam is divided into segments of equal length, the flexibility of 
each segment is the same, and the reduced form of the flexibility matrix of 
the elements (neglecting shear distortion) is given by Eq. 19. This may be 


4 10 0 0 0 0 
t @ 


compared with the unreduced form of the matrix shown in Eq. 12b. Note that 
the first and last row and column have been omitted because no moments are 
developed at the ends of the beam. 


Force Transformations: 
Moments in the beam over the supports at points 3 and 5 are taken as the 
redundants (Mx and Mx; ). On this basis, the force transformation matrices 


take the reduced form indicated in Eq. 20. 
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6 
mi 0000 0:6 o| |P 
\ 4 (20) 
0 3 0 013 3 P 
00000:0 6 P 
Mn 
0 00 3«3;0 3 X,| where 
é 1 
™, 0 0 0 6:10 O 
or 


\3 


[bei 


Results: 
Following the matrix operations of Eqs. A-1 through A-4, the results 
shown in Eqs. 21 were obtained. 


3 
DB= = (21a) 
1 
64 80 27 -9 12 
_ 228 | -16 -20 36 45 -60 (21b) 


784 296 -54 18 -24 
200 592 -108 36 -48 
-480 -162 54 -72 
b = | 7-144 -180 495-108 144 (21c) 
96 120 =216-270 360 
48 60 =108 549 -504 
0 0 01368 


15744 11472 =-2592 864 -1152 
12288 -3240 1080 -1440 


Er 4806 -1944 2592 (21d) 
(Symmetrical) 5778 -9072 
44928 


Eqs. 21c and 21d may be used to evaluate the forces and deflections in the 


beam produced by any system of applied loads R as shown by Eqs. A-6 and 
A-7. 


BS 
= 
r 
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IV. Continous Beam on Elastic Supports 
As an example of the procedure for treating flexible supports, the stress 


and flexibility matrices will be calculated for the beam system shown in 
Fig. 6. The beam is uniform throughout its length, and all supporting springs 


Fig.G Beam on Elastic Supports 


have the same stiffness. This example will also be used to demonstrate how 
a multiply redundant system may be analyzed by parts to avoid the inversion 


of a large matrix. The standard solution will be outlined first for compari- 
son. 


Standard Solution: 


For this uniform system, the flexibility matrices for the beam segments 
and the elastic supports are as shown in Eq. 22a. Thus the flexibility matrix 
of the complete system of elements takes the reduced form of Eq. 22b. 


fe 


fe 


Selecting the moments in the beam over the interior supports as the redun- 
dant forces, the force transformation matrices take the (reduced) form 
shown in Eq. 23. 

The flexibility and force transformation matrices of Eqs. 22b and 23 having 
been established, the system could now be analyzed by the matrix operations 
of Eqs. A. In the process, because this structure is indeterminate to the 
sixth degree, it would be necessary to invert a sixth order matrix D. This 
would offer no great difficulties in the present case, but for the illustrative 
purposes the analysis will be carried out here by parts. (It may be noted 
that considerable simplification could also be achieved in the analysis of 
this structure using the concepts of symmetry and anti-symmetry, but such 
ideas are not pertinent to the present discussion). 
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0 0000 0 
82000000000 0 0 
8 20000000 0 0 0} 
8 20000000 0 0} 
820000000 0 
06 0 8 
0 0 0 0 0 O (29) 
2 cet 3.00000 0| 
08 8 6 
3 000 0 
(Symmetrical) 3 0 
3 0 4 
3 0} 
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Analysis by Parts: 


For this example, the structure will be analyzed by dividing it into two 
parts, analyzing them separately, and then rejoining them. This symmetrical 
structure will be divided into two identical parts, as shown in Fig. 7; thus 
only one subsection need by analyzed. In the analysis it will be observed that 
the process of joining the two subsections together is identical with the 
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Fig. 7 Left-hand Subsection 


standard indeterminate analysis, except that the “primary” structure is 
indeterminate and therefore requires a preliminary analysis. 

The flexibility matrix for the subsection involves one element which is 
different from those representated in Eq. 22 b — the half length beam seg- 
ment at the right end. For this segment, the flexibility matrix is given by 
Eq. 24a. The complete flexibility matrix of the subsection elements then 


3 
fa (2) [F (24a) 


may be written in the reduced form of Eq. 24b. Selecting the beam moments 
over support 1 and 2 as redundants, the force transformation matrix takes 


0; 
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the reduced form shown in Eqs. 25. In Eq. 25a, the subscript “L” indicates 
that the moments and forces are referred to the left hand section. It should 
be noted that the forces P, and Mg which will be needed later to satisfy the 
continuity conditions between the two subsections have been included as part 
of the external force system acting on the subsection. Also it is seen that 
only two redundants are included in this subsection, hence its analysis re- 
quires the inversion of only a second order matrix. 
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The flexibility of the subsection, obtained by applying the matrix opera- 
tions of Eqs. A, is found to be as shown in Eq. 26. 


747 189 -27 -45 <-39 12 
459 243 -27 -153 -252 
459 189 -9 -396 

747 1065 636 
(Symmetrical) 2111 2236 
4064 


In order to join the subsections together, the flexibility matrix of the ele- 
ments of the complete structure is formed, as usual, as the diagonal matrix 
of the element flexibilities (subsections, in this case), as shown in Eq. 27. 
For this case, FR = Fy, because the two subsections are identical. 


© 
E | (27) 


The force transformation matrix required for combining the subsections 
takes the particularly simple form shown in Eq. 28 (this is typical of force 
transformations relating previously analyzed subsections). It should be 
noted how the force transformation takes account of the reversal of sequence 
in numbering between the left and right hand sections. In this equation, the 
external forces acting on the right hand subsection just to the right of the 
center joint are identified by “primes,” (e.g. A ), while the redundant inter- 
connection forces are designated Xp and Xm- 


Analysis of this system also requires the inversion of only a second order 
matrix. Carrying out the standard matrix operations on these flexibility and 
force transformation matrices, yields the flexibility coefficient matrix of the 
complete beam shown in Eq. 29. This matrix gives the vertical deflections of 


the beam due to unit loads, and is arranged in the sequence indicated by the 
external load matrix R of Eq. 28. 


™m, o'1 0 P 
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0:0 1 P, 

m: 110 0 

Pe 0 '-1 0 (38) 
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The stress matrix bg obtained (using Eq. A-3) during the joining of the 
subsections represents, of course, the relationship between the applied loads 


129.62 32.63 -4. 60 -6. 27 -1.87 0.12 
77-37 40.00 5444 =-3.28 =2.07 


F= 


76.33 38.59 4699 
74-41 38.00 4.99 

74-41 38.59 

metrical) 76. 33 


0.31 
-0.39 
-2.07 
-3.28 

5244 
40.00 
77.37 


0.06 
0.31 
0.12 
-1.87 
-6.27 
-4.60 
32.63 
129.62 


(28a) 


and the subsection force system, as shown in Eq. 30a. In order to obtain the 
forces in the beam segments of the subsections, it is necessary to go back to 


R 
- R 


(30a) 


the stress matrices by, and bp relating the subsection applied loads to the 


segment forces. This relationship is indicated in the matrix Eq. 30b. Thus 
the segment forces may be expressed in terms of the external load system 


(30b) 


applied to the entire beam as shown in Eq. 30c. Carrying out the operation, 
the stress matrix of the complete beam is obtained, as shown in Eq. 31. 
Only half of the matrix for this symmetrical structure is presented, in 


order to save space. 
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1359 -2175 307 418 125 
542 491 =-2053 473 468 15 
31 491 498 -2044 479 129 
310 487 <-782 -782 310 
8641 2175 =-307 -125 21 
2175 5158 2667 362 -218 -138 -26 
-307 2667 5089 2572 332 -223 -138 
-418 362 2572 4960 2533 332 -219 -125 


Naw 


where ™M, = 


DISCUSSION 


The matrix method presented herein provides a convenient and efficient 
means of calculating the flexibility and stresses of structures composed of 
beam type elements. After the basic flexibility and force transformation 
matrices of the elements of the structure have been set up (a relatively 
simple task for beam systems) the matrix manipulations required are of a 
routine nature which may be carried out by technicians having no training in 
structural analysis, or they may be handled in a routine manner by an elec- 
tronic digital computer. It should be emphasized, however, that an electronic 
computer is not at all an essential requirement of the method, however. All 
the calculations in the examples presented herein were carried out in 20 to 
30 hours on a desk calculator. 

One of the convenient features of this approach is that all the elastic 
characteristics of the structure are contained in the flexibility matrix, 
while the geometric characteristics are presented separately in the force 
transformation matrices. Thus changing the stiffness of any of the members 
requires only that the appropriate portion of the flexibility matrix be changed. 
Similarly, it is possible to determine the effects of shear and flexural dis- 
tortions separately, simply by using the proper flexibility matrices. Non- 
uniform members require the use of a special formula for evaluating the 
segment flexibilities, but introduce no other complications into the analysis. 

In the general description of this method given at the beginning of the 
paper, the statement was made that the inversion of the matrix D is 
equivalent to the solution of a set of simultaneous equations. Actually, of 
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course, the inversion of a matrix is equivalent to the solution of the equations 
for many different loading conditions. Consequently, if only a single loading 
condition is contemplated for the beam system usder consideration, it will be 
simpler to solve the simultaneous equations than to obtain the inverse of the 
matrix. In such a case, the inversion operation indicated in Eq. A-2 may be 
considered merely as a symbol of the solution of the set of equations. On the 
other hand, if the complete set of flexibility coefficients is required (as for 
example, for subsequent use in a vibration analysis) then the actual inversion 
operation must be performed. Langefors has presented a thorough discussion 
of this point in Reference 9. 
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APPENDIX 


Note on the Displacement Matrix Equations 


The dual function performed by the force transformation matrix bj in 
Eq. A-1 (or bp in Eq. 13c) may appear somewhat surprising at first. How- 
ever, the validity of the expression may be demonstrated easily by consider - 
ing it as a statement of the principle of virtual work. By this principle, the 
deflection at a specified point of a structure, may be obtained by applying a 
unit load at the location and in the direction of the desired deflection and ob- 
taining the product of the resulting force in each member multiplied by its 
specified deformation. The sum of these products for all members in the 
structure then equals the desired deflection. 

In Eq. A-1, the matrix product Fb; represents the internal deformations 


of the elements, while the matrix b} is the (transposed) matrix of forces in 


the members due to unit loads applied individually at points where the deflec- 
tions are to be found. The matrix product b) (fb ,) thus gives the sum of the 


member force-deformation products as required to obtain the desired deflec- 
tion. The matrix by enters the expression twice because deflections are 
being calculated at the same points as the loads are being applied to. 

It may also be worthwhile to note that Eq. A-1 is a typical equation of co- 
ordinate transformation expressed in matrix form. (See for example 
“Elementary Matrices” by R. A. Frazer, W. J. Duncan, and A. R. Collar, 
Cambridge University Press, 1955, p. 29). In this case, the transformation 
is from the deformation coordinates of the individual elements as exempli- 


fied by f to the deformation coordinates of the combined (primary) system 
represented by D. 
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PLASTIC DESIGN OF COVER PLATED CONTINUOUS BEAMS 


E. P. Popov* A.M. ASCE and J. A. Willis** J. M. ASCE 
(Proc. Paper 1495) 


SYNOPSIS 


A non-prismatic continuous beam may be considered well designed by the 
plastic method if at the ultimate load a plastic hinge forms throughout the 
beam. Significantly, according to the plastic method of design, an infinite 
number of alternative designs to carry the same applied load is possible. 
The material of a beam may be distributed quite differently to achieve the 
same ultimate carrying capacity. This situation is entirely different from 
that found in elastic design where theoretically, once the type of cross- 
section is fixed, there is only one solution of a given problem. In this paper 
alternative plastic designs for a two-span continuous beam loaded with con- 
centrated forces are discussed. Experimental evidence is submitted to serve 
as a justification for the simple plastic method of design of non-prismatic 
members and to show some of the shortcomings in the simple plastic theory. 


INTRODUCTION 


For ductile materials, such as structural steel, plastic methods of analy- 
sis and design are applicable. Realistic estimates of the carrying capacities 
of structures may be obtained on this basis. Convenient analytical procedures 
have been developed(1,2) for the design of structures based on this method. 
Reliable experimental verification of the theory is also available(3) for many 
practically significant cases. However, so far some aspects of the problem 
received insufficient emphasis. One of these deals with the fact that the plas- 
tic method of design permits an infinite flexibility in the manner of propor- 
tioning the various members of a structure. The material may be distributed 
differently to achieve the same ultimate carrying capacity. For example, ina 
simple portal frame fixed at the base a variety of column and beam sizes may 
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be shown to be adequate for a given loading condition. By making columns 
stronger the beam may be reduced in size, or vice versa, provided of course, 
that the combined plastic strength of the system suffices to carry the applied 
load. Several adequate practical designs are possible even though the beam 
and columns are limited to constant cross-section throughout their respective 
lengths. An optimum design to achieve the minimum weight of such a portal 
frame has been discussed. (4) 

Greater possibilities open up for the weight optimization of a structure if 
the cross-section of each member is permitted to vary along the length. For 
I-beams and wide-flange sections this is effectively accomplished by means 
of cover plates. Theoretical possibilities of this variation for uniformly 
loaded fixed-ended beams has been considered.(5) In the interest of weight 
economy the practical necessity for varying the cross-section of beams and 
columns is apparent in the recommended design procedures.(6,7,8) This pos- 
sibility of having numerous designs for indeterminate structures to choose 
from when the plastic design method is used is unique. According to the elas- 
tic theory, whether the structure is determinate or indeterminate, there is 
one and only one solution of the problem for a given load system. 

In this paper, several designs of a two-span continuous beam loaded by 
concentrated forces in the middle of each span are considered first. A some- 
what analogous case of a three-span continuous beam is also discussed. Then, 
since most of the experimental research in plastic design of continuous beams 
and frames has been confined to prismatic members, experiments on the de- 
signed beams were performed. These experiments may serve in some mea- 
sure in justification of predictions of the ultimate capacity of non-prismatic 
members based on the simple plastic theory. 


Design of a Continuous Beam 


As an example of plastic design of a continuous beam, consider the beam 
loaded as shown in Fig. la. Here it is evident that if the beam were not con- 
tinuous over the center support, the moment diagram in each span would be 
triangular in shape with a maximum ordinate of +PL/4. The same diagrams 
may be considered as a part of a total bending moment diagram for a con- 
tinuous beam, except that for a continuous beam the base line is a line ABC, 
Fig. 1b. 

It should be noted particularly that the ordinates PL/4 in each span are 
based entirely on the conditions of statics and that irrespective of the mag- 
nitude of P these expressions remain the same. In plastic design, the ordi- 
nate for the moment diagram at the center support may be chosen arbitrarily. 
In terms of the notation shown in Fig. lb, this means that the value of a di- 
mensionless parameter o@ may be chosen arbitrarily. Actually, @ may vary 
between the values of zero and two.(9) At both limits the continuous beam 
ceases to exist. For @ = 0, two simple spans are obtained; for a@ = 2, the 
beam becomes a cantilever. 

After a definite value of a is selected, the moment diagram becomes 
fixed. If, further, P is considered as an ultimate load, P,,)¢, sections along 
the beam may be proportioned following the rules of plastic design for the 
corresponding values of My. A beam of continuously varying size selected 
on this basis would be efficient from the point of view of plastic design. At 
an ultimate load a hinge continuous across the entire span would form all at 
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once. In an actual design of such a beam, portions of the span where the 
bending moment is small must be modified to care for the shear require- 
ments. 

The arbitrary choice of @ is permissible only on the basis of plastic an- 
alysis. For the beam considered, when material behaves elastically, the 
general shape of the moment diagram remains the same. However, there is 
only one particular value of a@ which will satisfy the conditions of continuity 
at the middle support. The proportions of the beam influence the magnitude 
of a. A beam made stiff over the center support will develop there a large 
moment; a beam which is flexible over this support would develop there a 
smaller moment. Theoretically it is possible to select a continuous beam of 
constant strength for a given shape of a cross-section. This solution would 
be unique.(10) By contrast, for the same conditions an infinite number of de- 
signs is possible according to the plastic theory. Therein lies one of the 
fundamental differences between the elastic and plastic methods of design. 

In structural work where rolled sections are normally used, a continuous 
variation in the size of a beam is not practical. Instead, a certain minimum 
size beam may be selected and by using cover-plates the bending moment re- 
quirements may be approximated step-wise. On this basis four two-span 
continuous beams of different cross-sections to carry the same loads P were 
designed and tested. 

One of the continuous beams was made prismatic. For such a case the 
moment resisting capacity of the beam at the applied forces P and at the cen- 


ter support is the same, i.e. a = :. Hence equating moment at the applied 


load to the moment at the support, see Fig. 1b, one finds 


Pint 6M 
- 1/2 M, = M, or P 
Pp u 


_ | 
where 
1S the ultimate transverse load 
M, is the ultimate plastic moment 
L is the length of each span 
Three additional beams were made with cover-plates in which the plastic 


moments were different at various sections. For the loading conditions con- 
sidered, the basic equation, a generalization of Eq. 1, is 


Pat” 
4 = Moa + 1/2 Mab 


where Mpa is the plastic moment at the applied loads 
Mpb is the plastic moment at the middle support 


A 6 in. 12.5 lb. I-beam 16 ft. long beam was selected as the prismatic 
specimen. Based on the assumption that the yield strength of structural steel 
is 40 ksi, from Eq. 1 it was found that Pyjt should be 20.9 kips and that P,)j¢ 
L/4 is very nearly 500 in. kips. 

A 5-in. 10 lb. and a 5 in. 7.7 lb. continuous I-beam were coverplated in the 
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middle of each span to provide for very nearly the same P,,); as was deter- 
mined for the 6 in. prismatic beam. With no cover plates over the middle 
support, Mpb is readily computed, and, from Eq. 2, the Mpa; the bending mo- 
ment capacity at the applied loads, was found. By assuming the same yield 
strength for cover plates as for the rolled sections, the appropriate sizes of 
cover plates were determined. Another 5 in. 10 lb. I-beam with cover plates 
only over the middle support was designed to carry the same Pyjt. The de- 
sign results are summarized in Table I. 


TABLE I 

Member Cover Plates (kip-in) (kip-ir.) P 

Size Size (in.x.in.) plates plates ult 
cips in (kips ) 
kips- 

61 12.5 none 334 --- 501 20.9 

§110.0 3/l6éx4 222 378 489 20.4 

41 7.7 3% and} x3 138 430 499 20.8 

5110.0 7/16 x 3% 222 555 500 20.9 


The beams selected were so chosen that shear would not be important. 
Further, assuming that 1.88 is the proper factor of safety(11) for obtaining 
working loads from the ultimate loads, it was found that at working loads in 
all cases the deflections of the selected beams did not exceed 1/360-th of the 
span. However, for the 4 in. 7.7 lb. beam having cover plates in the middle 
of each span (Beam No. 5) yielding over the center support takes place before 
the total working load is reached. In plastic designs such situations are 
possible. 

A composite bending moment diagram for the four cases considered is 
shown in Fig. 2. A different baseline is applicable in each case. From such 
a diagram cut-off points for cover plates were determined. As a general 
rule, the cover plates were extended 6 in. beyond the theoretically required 
cut off points. 

Besides the four continuous beams described above, another related con- 
tinuous beam was analyzed and tested. The loading diagram for this three- 
span continuous beam is shown in Fig. 3a. If the end spans of this beam were 
rotated downward at B and C through 90°, a small portal frame would result. 
Hence this beam simulates approximately an unfolded portal frame. On the 
other hand, if the loading diagram is viewed upside-down, the loading arrange- 
ment resembles the case considered earlier. The principal difference be- 
tween the two cases is in the shape of the collapse mechanism, compare 
Fig. lc and Fig. 3c. For this case Eq. 2 does not apply, instead, as may be 
readily seen from Fig. 3b the basic equation should be 


_+M (3) 


Where Mpa is the plastic moment at the applied load 


Mpb is the plastic moment at the interior supports. 
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A 5 in. 10 lb. I-beam with 7/16 in. by 3 1/2 in. cover plates which exactly 
duplicated on the previously designed beams having Mpa = 555 kip-in. and 
Mp} = 222 kip-in. (See Table I) was selected for the experiment. For this 
case the cover plates extended about 1/2 in. beyond the required cut off points 
based on plastic analysis. Significantly, the elastic analysis indicated first 
yield would take place at the edges of cover plates. Therefore, the elastic 
and plastic analyses lead to different requirements for the length of cover 
plates. Further research on this question appears to be necessary. 

For control purposes, a 5 in. 10 lb. I-beam to span 3 ft. 6 in. and a 5 in. 

10 lb. I-beam with 3/16 in. by 4 in. cover plates to span 4 ft. 2 in. were se- 
lected. These beams were made short in order to obtain gradients in the 
bending moments which would be comparable to those found in the selected 
continuous beams. 

Details for the seven beams chosen for the experiments are shown in 
Fig. 4. Note that 5/16 in. stiffeners were provided at all supports and at all 
points of concentrated load application. A liberal amount of intermittent weld- 
ing was specified to provide interconnection of all beam elements. 


Experimental Procedure 


All beams were fabricated by a local steel fabricator and were made of 
ordinary structural steel. Beams of the same size were made from the same 
lot and all coverplates of similar size were likewise taken from same lots. 
Beams were not annealed either before nor after welding. 

After the specimens were delivered to the laboratory, cross-sectional di- 
mensions were measured. At a later date, unreinforced portions of the beams 
were cut out and weighed. No significant discrepancy with the values given in 
the AISC Manual were found. Therefore, handbook .values of section proper- 
ties were used in the calculations. 

To verify the assumed mechanical properties of the material, 1 1/2 in. by 
8 in. gage length tensile coupons were made. Three coupons were taken from 
each size I-beam (one from each flange and one from the web) and one from 
each size of coverplate. 

A study of the stress-strain diagrams for the coupons showed no anomali- 
ties, although well defined yield strengths varied from 32.4 ksi to 45.4 ksi. 

An average value for the yield point for a given cross-section was determined 
by using a weighted average based on the area and location in the beam. The 
average values used in the final analysis are tabulated in Table IL. 

The beams were tested in a 200,000 lb. Olsen testing machine which has a 
long horizontal bed plate as may be seen in Fig. 5. The two experimental 
simple beams were loaded at midspan. The four two-span continuous beams 
were loaded symmetrically at midspans with a concentrated load. This was 
accomplished by applying, by means of the machine of a single concentrated 
load to a heavy distributing beam. A spherical block was used between the 
distributing beam and the head of the testing machine. The arrangement for 
the three-span experiment is shown in Fig. 5. The central span was loaded by 
the machine in the middle; the end supports were not permitted to deflect by 
maintaining pressure in calibrated hydraulic jacks, (see Fig. 6). A record of 
the required pressures in the jacks provided information on the magnitude of 
end reactions. 


In all experiments, beams were placed on rocker supports. Spherical 
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blocks and rollers were extensively employed. In addition to plates and 
shims, Hydrostone cement was applied at joints to assure good bearing sur- 
faces. 

To prevent sidewise buckling of the specimens at midspans of beams, long 
wire ropes were attached to the vertical bars which in turn were attached to 
the beam stiffeners, as may be seen in Figs. 7 and 8. 

Deflections were measured at midspan in all cases except in the end spans 
of the three-span beam. In most instances such deflections were measured 
by means of dial gages as may be seen in Fig. 8; in others when a Sanborn 
Recorder was available, deflections were obtained by Bourns Linear Poten- 
tiometer. This set-up is shown in Fig. 7. For the two-span beam experi- 
ments some discrepancy between the deflection readings in the two spans 
were observed. The two readings varied from 3% to 22%. However only the 
average values of deflections are reported since it is believed that this dis- 
crepancy does not affect substantially the end results. 

In several experiments strain measurements were made on the flanges on 
the side opposite the point of load application. These were recorded with the 
aid of the Sanborn Machine. In all cases the initial yielding occurred prior 
to the predicted values. It is believed that residual stresses may account for 
this observation.(12,13,14) The data obtained are not reported. 

The appearance of the beams at the end* of the experiments may be seen 
in Fig. 9. The experiments usually were discontinued after the deflection 
exceeded 1/36-th of the span length. The wire rope bracing was considered 
unsafe beyond such a point. In any event, since a substantial leveling of the 
load-deflection curve occurs in this range of deflection, continuation of the 
experiments beyond this point, without much additional load, would precipi- 
tate some local collapse. 

Prior to the load application, the beams were white-washed. This aided in 
the detection of yield lines which occurred during loading, (see Fig. 9). 


Experimental Results and Conclusions 


In Fig. 10 are plotted load-deflection diagrams for beams 1, 2, and 7; in 
Fig. 11, for beams 3, 4, 5, and 6, each of which according to the plastic 
theory should carry practically the same ultimate load. The load at the time 
when maximum deflection reached 1/36-th of the span length was used to de- 
fine the ultimate load, Pit: This amount of deflection is ten times greater 
than is currently permitted by AISC specifications. With the exception of 
beam 2, where the strength appears to decay around this amount of deflection, 
all beams still could have carried slightly higher loads. However, the lateral 
bracing was considered inadequate to continue loading until actual collapse of 
the beams would occur. 

On the graphs in Figs. 10 and 11, in the plastic range, two points are usu- 
ally shown corresponding to the same deflection. The higher points are for 
initial reading of the applied load; the lower points indicate the applied load 
after a lapse of time, which permitted the material to complete plastic flow. 
All reported ultimate loads P.it are based on inter- and extrapolation of the 
curves passing through the lower points. 

The principal results of the experiment as well as some derived quantities 
are summarized in Table II. For further clarification, the significant results 


* These beams were tested and later re-tested, to obtain the load- 


deflection data reported in this paper. 


= 
oe 
is, 
uw 


ASCE PLASTIC DESIGN 1495-7 


are plotted on the bar graph of Fig. 12. The load carrying capacity of beams 
at the actual yield point strength of steel has been computed by conventional 
methods. The yield strength of steel was based on data obtained from tests of 
coupons. 

Based on the results of these experiments and the analysis of the data, a 
number of conclusions appear to be warranted: 


a) As anticipated by theoretical considerations, it is possible to design 
non-prismatic beams using the simple plastic theory. 


b) Although the simple plastic theory appears to have led to a very con- 
sistent design of beams 3, 4, 5, and 6, due to strain hardening of the material, 
the real capacity of beams is higher than such a theory predicts. In litera- 
ture dealing with plastic analysis it became customary to cite cases showing 
remarkable agreement between experimental results and predictions based 
on simple plastic theory. Such one-sided exposition does not seem to be en- 
tirely proper. The simple plastic theory expands our understanding of struc- 
tural action and leads, in a very simple manner, to reasonable estimates of 
the ultimate capacity of structures. It should be recognized, however, that 
this theory does not give the complete picture. 


c) In examining Fig. 12 it should be noted that the elastic theory is com- 
pletely unable to predict the behavior of beams beyond the yield point load. 
The regions cross-hatched in the diagram simply cannot be even estimated 
using the elastic theory. A substantially equal ultimate performance of beams 
3, 4, 5, and 6, may be contrasted with a rather varied predicted performance 
of the same beams based on elastic analysis. Plastic design should, there- 
fore, promote greater consistency in design for strength. 


d) It is significant to reiterate that in plastic design for the same ultimate 
load capacity different size members may be utilized. Material may be dis- 
tributed quite differently to achieve the same ultimate carrying capacity. This 
has practical implications of economy by permitting a freer choice of mem- 
bers. Also, the total weight of construction may be minimized. 


e) Direct experience with testing beams into the plastic range forcefully 
pointed to the dire necessity of bracing beams laterally. Therefore, it seems 
safe to say that in actual design, likewise, to assure the desired performance 
of a structure, at locations where plastic hinges may be anticipated, bracing 
must be introduced. Such a general rule seems wise to apply whether one 
designs using elastic or plastic method. 


{) In the reported experiments, coverplates performed satisfactorily. 
However, further research is needed to determine the minimum amount of 
welding and the proper length of the coverplates for plastic design. 
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(a) Beam Loading Diagram 


(6) Moment Diagram 


(c) Collapse Mechanism 


Fig. 
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2 
(2) Beam Loading Diagram 


Mpa 


(6) Moment Diagram 
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(c) Collapse Mechanism 
Fig. 3 
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SEA BOTTOM PRESSURE FIELDS PRODUCED BY YAWED VESSELS 


P. M. Fitzpatrick* 
(Proc. Paper 1496) 


ABSTRACT 


An understanding of the flow pattern about an ellipsoid immersed at rest 
in a parallel stream is often of value in aerodynamic and hydrodynamic 
studies. The steady flow, without circulation, about an ellipsoid immersed at 
rest in an infinite half-space of perfect fluid is discussed. A uniform stream 
flow is considered in which the centroid of the ellipsoid remains at a finite 
distance h from the plane boundary of the half-space. The major axis of the 
ellipsoid is taken to be parallel to the plane boundary and to make an angle @ 
with the direction of the undisturbed flow. Approximation equations are given 
for the velocity and pressure distributions on the plane boundary. These 
equations are based upon the method of images and the well-known theory of 
the flow pattern about a stationary ellipsoid in an infinite fluid medium. A 
procedure for obtaining the lines of equipressure change on the boundary is 
also described. The sea bottom pressure distribution produced by slowly 
moving yawed surface vessels in advance of the upstream hull-line may be 
approximated from computations based on formulae developed for the ellip- 
soid model. These approximations should hold when the ship proceeds 
through calm water of moderate depth at an angle of yaw sufficiently large to 
produce a wake, say 6 greater than 20°. As an illustration of the application 
of the theory certain results of calculations of sea bottom pressure distribu- 
tions and lines of equipressure change are given for a representative case. 


INTRODUCTION 


In certain areas of hydrodynamics it is important to be able to predict 
quantitatively the sea bottom pressure distribution associated with a surface 
ship or submarine moving slowly** in a calm sea. The most frequently 


Note: Discussion open until June 1, 1958. A postponement of this closing date can be 
obtained by writing to the ASCE Manager of Technical Publications. Paper 1496 is 
part of the copyrighted Journal of the Engineering Mechanics Division, Proceedings 
of the American Society of Civil Engineers, Vol. 84, No. EM 1, January, 1958. 

* Hydrodynamics Branch, Research Div., U. S. Navy Mine Defense Lab., 
Panama City, Fla. 

** A conservative criterion for the validity of an assumption made later that 
the surface behaves as a rigid reflector is that the vessel proceed at a 
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encountered case is, of course, the one in which the longitudinal axis of the 
vessel moves approximately parallel to the track (i.e., the angle of yaw, 6, 

is less than say 3°). For this case, which is amply treated in the literature, 
the pressure distribution functions are usually obtained either from a study 
in which the ship is represented by the source-sink technique or from a study 
of the flow pattern past an ellipsoid whose dimensions have been chosen to 
simulate the actual ship profile. Little information appears in the literature, 
however, on the methods which one may use to approximate the sea bed pres- 
sure distributions when a slowly moving surface (or subsurface) vessel pro- 
ceeds at a constant angle of yaw greater than 3°. These latter cases, although 
encountered infrequently, do on occasion find practical application and hence 
merit attention. 

When a vessel proceeds through a calm sea at low speed, the magnitude of 
the angle of yaw, 0, affects the type of flow pattern which occurs. For values 
of 6 less than approximately 3° the usual treatments of the problem assume 
that a pure potential flow without circulation will hold. For values of 9, in the 
range 3°< 6 <15° to 20°, an approximation to the actual flow pattern about 
the ship profile may be obtained if one supplements a pure potential flow by 
the inclusion of circulation. This type of flow pattern is not treated in this 
paper. For values of 6 which are greater than 20°, a wake forms aft of the 
vessel and the methods of analysis presently available do not permit one to 
make quantitative predictions of the pressure distribution in this aft region. 
Nonetheless, a pressure distribution function obtained from pure potential 
flow theory (without circulation) should hold, as it does in the case that @ is 
less than 3°. The region of validity of solutions based on this assumption 
must, however, be restricted to the neighborhood of the upstream hull line. 

If 6 is greater than 20° it is possible to obtain a sea bottom pressure dis- 
tribution function, which should hold approximately in the forward region, 
from a consideration of the pure potential flow (without circulation) about an 
ellipsoid. To this end a treatment of a yawed ellipsoidal body immersed at 
rest in an infinite fluid half-space has been made. Expressions which may be 
used to compute the velocity and pressure distributions on the plane boundary 
of the half-space are given. In addition, the problem of determining the 
curves of equipressure change on this boundary has been considered and a 
possible method of constructing these curves is developed and discussed. 


Theoretical Procedures 


Consider an ellipsoid form with semi-principal axes a, b, and c, (a>b>c). 
Let it be immersed at rest in an incompressible, nonviscous, infinite fluid 


medium which is moving with uniform total velocity y uv? + v2. An irrotation- 
al flow is assumed. Let the reference frame be fixed in the ellipsoid and so 


oriented that the Z-axis is perpendicular to the velocity vector v2 + v? , 
The major axis of the ellipsoid (i.e., the X-axis) is assumed to make an angle 
6 with the direction of the uniform stream. Unless stated otherwise, in what 
follows, it is assumed that @ satisfies the condition 90° > @ > 20°, The 
velocity components of the uniform stream are designated as -U and -V 


speed (¥ gh, where h is the depth of the water and g is the magnitude of the 
acceleration of gravity. 


at 
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respectively and are directed along the negative X and Y-axes as shown in 
Figure (1). 
; For the system just described, the velocity potential ¢ , is given as(1) 


y 
UX+ VY; abcUX 2 abcVY dA a) 
2—a, (a“+A)ky 2-B, (b“+A)k y 


Here A is an ellipsoidal coordinate which is related to the rectangular 
coordinates X, Y, and Z through the equation 


9 9 a 2 
(X/a) J (Y/a) (Z/a) (2) 
1+A/a2 b2/a24A/a2 


The symbols kj, @, and 8, which appear in Equation (1) are defined as 
follows: 


ky = a? (144/02) (b2/a24A/a2) ark’, 


(3a) 
(b/c) (c/a) 
c c/a 
(14A/a2) k’, (3b) 
d(A/a2) 
B,, = (b/a) (c/a) (3c) 


The elliptic integrals which appear on the right-hand side of Equation (1) have 
been tabulated. (2) 

Over the range of definition, Equation (2), regarded as a function of A, 
represents a family of real confocal central ellipsoids. A nest of expanding 
quadric surfaces may be generated by successive substitution in Equation (2) 
of a monotonically increasing sequence of values of A. In particular, for A = 
0, the quadric surface coincides with the elliptic form shown in Figure (1). 

If the velocity components in the X, Y, and Z directions of the disturbed 
flow are designated by u, v, and w respectively, one finds by differentiation 
of Equation (1) that 


(X/a) K)(X/a) 


u = —06/dX = —U + 2 (b/a) (c/a) 
(1+A/ A, +A/a*) 


Ky (Y/a) d(r/e”) 
— (b/a) (c/a) Ky 
(b?/a2+A/a”) 9 


(4a) 
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(¥/a) 
v = —dp/dy = — V + 2 (b/a) (c/a) 


(14A/a2) 


(4b) 
Ky (Y/a) A/e”) 
(b2/a2+A/a2) / 2 (b A 
a 
(Z/a) K, X/a Ky (Y/a) 
w = —0$/0Z = 2 (b/a) (c/a) + (4c) 
The symbols Kj, Kg, and Ap are defined by the equations 
U 
K,= (5a) 
2-a, 
Ky = 
2-B, (5b) 
(X/a)” (¥/a)” (Z/a)” 
= — 
2 
(5c) 


Equations (4a), (4b), and (4c), which are valid only for an infinite medium, 
may not be applied directly to the problem of an ellipsoid in an infinite fluid 
half-space. They must be modified to take account of the fact that the normal 
component of the fluid velocity across the plane boundary must vanish. This 
latter requirement may be satisfied if one considers the system to be com- 
posed of an ellipsoid and its image in the boundary. The introduction of the 
image while it satisfies the velocity criterion on the plane boundary results 
in a perturbation of the flow at the ellipsoid wall and on the plane Z = 0. If 
the ellipsoid is not too close to the plane boundary, however, this interaction 
effect is small. For the purpose of a rough approximation, this perturbation 
effect may be neglected* so that one may obtain the velocity function on the 
plane boundary by the simple process of doubling the velocity changes given 
through Equations (4a) and (4b). 

If the plane Z = +h, h> cc, is taken to be the boundary, the modified expres- 
sions for the velocity components on this plane become 


* If a higher order of approximation is desired, a greater number of images 
may be employed. The general scheme for the positioning of images 
whose centroids lie on the Z-axis is given by the expression Z = t 2nh 


; 
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(X/a) K(X a) a 
Uy, —U +4 (b/a) (c/a) — 
(1+A Aj, (1+A/a*) b~/a~+A/a~) (6a) 
dA a’) 
—2 (b/a) (c a) Ky ———— 
>» 
(Y/a) a) Ko(¥/a) 
(b“/a“+A/a~)k (1+A/a*) (b2 asa a’) 
(6b) 
dla a2) 
—2 (b/a) (c/a) Ko neaimmnanenn — 
= 
» 
w, = 0 (6c) 
where A,= + + —— and the additional 


subscript h on the velocity components denotes that they are to be evaluated 
on the plane Z = +h. 

If the ellipsoid is immersed in a body of water sufficiently deep so that 
the effect of the surface may be neglected, Equations (6a), (6b), and (6c) may 
be applied to compute the velocity field on the sea bed forward of the ellip- 
soid. Equations (6a), (6b), and (6c) will also serve to describe the velocity 
distribution produced on the bottom in the forward region by a half ellipsoid 
moving in the surface of a body of calm water of moderate depth h if the 
plane Z = 0 may be treated as a symmetry plane. This latter condition is 
met approximately if the half ellipsoid is moving slowly through the fluid 
since for this case the surface of the water is essentially smooth. In many 
cases of practical interest which involve surface vessels it is possible to 
make a rough representation of the underwater hull by a half ellipsoid of 
suitable dimensions. When such a representation is used and within the 
limitations stated, Equations (6a), (6b), and (6c) may be used to compute the 
approximate sea bottom velocity distribution in the neighborhood of the up- 
stream hull-line. 

If the parameter a is chosen as the unit of length, one may for fixed values 
of the ratios a/c, b/c, and h/a express Equations (6a), (6b), and (6c) ina 
functional form which will prove to be convenient for the purpose of the sub- 


sequent discussion. If note is taken of the fact that V = U tan 6, one may 
write 


uj, U (X,Y,A) (7a) 


v, = U [B (X,Y,A) -tan 6] (7b) 


where a@ and B are defined by the relations: 
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X xX 
a = a(X,Y,A) = 4 (b/a) (c/a) 
(1+A/a“)k', Ay, (2-a,) (1+A/a“~) 
(8a) 
(Y/a) tan 6 2 d (A/a2) 
(2-8) (b“/a~+A/a“~) (2—a,) 
(Y/a) (X/a) 
B = B(X,Y Abed (b/a) (c/a) 
Ay (2—a,) (1+A/a") 
(8b) 
Y/a) tan 0 2 tan d (A/a?) 
(2~8,) (b“/a“+A/a*) (2-B,) 
/a“ 
Although A is related to X and Y through the relation 
(X/a) (Y/a) (h/a) (9) 


(14A/a2)— (b2/a24A/a”) (c2/a24A/a2) 


it is convenient to think of a and # in the functional form given in Equations 
(8a) and (8b). Here, and in what follows, the parameters U and @ are con- 
sidered to be fixed and, unless stated expressly otherwise, the velocity and 
pressure fields will be referenced to the plane boundary. Thus the velocity 
field at any arbitrary point (X', Y', h) is determined by the triplet of values 
(X', Y',A') which satisfies Equation (9). 

The pressure distribution function Ap, may be obtained from the Bernoulli 


equation for steady flow. For the case under consideration, the Bernoulli 
equation becomes 


Ap = pU2/2 [1 + (10) 


In Equation (10) the units of Ap will be pounds per square foot when the 
velocities are expressed in ft/sec and the mass density, P, is expressed in 


slugs per cubic foot. It is more convenient to express Equation (10) in the 
dimensionless form, 


AP = tan 9) (11) 


pU2/2 


in which, the squares of a and § have been neglected because a and § are 
small compared to unity. 

In principle, a point-wise description of the pressure field has been given 
since the pressure change at any point (X, Y, h) may be calculated from 
Equation (11) after a and 8 have been determined from Equations (8a) and 
(8b). Of equal practical interest, however, is the determination of the lines 
of equipressure change. In the discussion which follows a procedure for 
obtaining these lines of equipressure change is described. 
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It is clear from an inspection of Equation (11) that at each point which lies 
on a line of given equipressure change, the following condition must be 
Satisfied 


tan 8) = K 


(12) 


where K is a constant which has the same value at each point on the line. For 
a prescribed value of 9, one may write Equation (12) in functional form as 


f (X,Y,A) = K (13) 
If one were to express A as an explicit function of X and Y, say A = g(X,Y), 
by solving Equation (9) for A, and then were to substitute g(X,Y) into Equation 


(13), the following relation would be obtained 


f (X,Y,g (X,¥)] = K = F (X,Y) (14) 


where F is some other function of X and Y only. One might expect that 
F(X,Y) would represent a complicated curve in the Z = +h plane. Nonetheless, 
it would appear that the equation of a line of equipressure change could be 
easily constructed from the set of coordinate pairs obtained when one sub- 
stitutes arbitrary values of Y into Equation (14) and solves for the corres- 
ponding values of X. While the foregoing development may be carried out in 
principle, a careful study of Equations (7a) and (7b) will show the difficulties 
which are inherent if one proceeds in such a fashion. These difficulties arise 
because of the complicated form of the expressions for a and § and because, 
for points which do not lie on the coordinate axes, Equation (9) represents a 
cubic equation in A. 

On the other hand, it is possible to express Equation (12) as a quadratic 
equation in Y which has as its solutions 


(15) 


where 1, s, r, and t (to be defined) are functions only of A and K, or A alone. 
In terms of Y, Equation (9) yields 


Y=: [m—nX2] (16) 


where m and n are functions only of A to be defined explicitly below. If Y is 
eliminated between Equations (15) and (16) one has, after performing two 
squaring operations to remove the radicals, a quadratic equation in x2 


(21 (mnl—ms+nt) + (mr“—2st) } 


(nl2422—2s!) + (17) 


21 (mnl—msent) + (mr2—2st) —4 (mist)? + 87] 
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The functions 1, s, r, t, m, and n introduced above are defined as follows: 


1 6 d 6 
le 


-1/8 le >) k', K 
1/8 (a/b) (a/c) ky (18a) 


9 


to 
| 
t 
| 


tan 6 
+ 

9 9 
a*(1+A/a) (b2/a7+A/a”) 


P 
s = - + — 1/8 (a/b) (a/c) k') K (18c) 


a-(14A/a2) 


9 
(h/a)“ k’ d 
+ 1/8 (a/b) (a/c) —- {1 8d) 


9.,9,9 9,2 2 ( 2 2 A/ 2) 
/ 9 a~ c“/a"+A/a 
é 


(h/a)” 
(c2 az+A/a2) (18f) 


maa (b2/a2+A/a2) 


The symbols d and e which appear in Equations (18a), (18c), and (18d) are de- 
fined by the relations: 


d (A/a?) 


2 (14A/a7)k', 


d (A/a”) 


(b2/a74A 


(19b) 


For a chosen value of A and of K there are four values of X which satisfy 
Equation (17). These four values of X when substituted into Equation (16) 
serve to determine four corresponding values of Y. Thus to a single value of 
A and of K there correspond four points of equal pressure change. By making 
successive choices of A for a fixed value of K one can obtain a sufficient num- 
ber of points of equal pressure change to use as the basis for the construction 
of a given line of equipressure change. If this procedure is repeated one may 
obtain, by varying K, any desired number of lines of equipressure change. 
Certain limits must be placed on the value of A which may be chosen if the 

four values of X obtained from Equation (17) are to be real. This restriction 


1 
(18b) 
sl 
2-a, 2 2-a,, 2-B, 
i 
Tt 
b2 24A/a2) 
“/a~+A/a™ 
n= (18e) Ve 
(14A/a*) 
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on the permissible values of A arises because a given curve of equipressure 
change will not, in general, intersect all of the ellipsoidal ecm which 
would be generated if A were permitted to vary over the range, -c 2< <x. 
If the line of equipressure change neither pierces nor is tangent to the ellip- 
soidal surface which corresponds to a chosen value of A, the values of X ob- 
tained from Equation (17) will be imaginary. 

For values of A less than (h2 - c2), the confocal central ellipsoidal sur- 
faces defined by Equation (2) will not intersect the plane Z = +h. Hence, for 
such values of A, the roots of Equation (17) will be imaginary. When A satis- 
fies the condition that A is greater than (h2 - c2), the confocal central ellip- 
soids intersect the plane Z = +h in the real closed elliptic curves given by 
Equation (9). If now the value of A is so chosen that the given line of equi- 
pressure change intersects the real closed elliptic curve associated with the 
value of A, the roots of Equation (17) will be real. 

The findings from theoretical studies and experimental measurements of 
flow patterns around ships will serve as a guide to the general location of the 
positive and negative regions of the pressure fields associated with ships. 
Thus, for a particular case under study, if an ellipsoidal representation of the 
vessel is used, it is a simple matter to choose the proper range of values of 
A so that the roots of Equation (17) are real. 

The results of certain computations which have been based on the use of 
the ellipsoid model, a:b:c = 10.50:1.11:1, are illustrated in Figure (2). The 
ratio h/a was taken to be 0.35. The ellipse shown in Figure (2) represents 
the projection on the plane of the bottom of the principal horizontal cross- 
section of the ellipsoid model. 

Computations have been carried out for three different angles of yaw, (6 = 
30°, 45°, and 60°) and the line of zero pressure change to be associated with 
each angle is shown on Figure (2). For the value 6 = 45° the dimensionless 
quantity 4 p/p (u2 + V2)/2 has been evaluated at a number of points which lie 
on the bottom in the region where the theory may be expected to apply. At 
each such point the value of Ap/? (U2 + V2)2 has been set forth on the figure. 

Equation (9) defines the curves of intersection of the confocal central ellip- 
soidal surfaces given by Equation (2) with the plane Z = +h. These curves are 
real central ellipses (with respect to (0, 0, +h)) if A satisfies the condition 
that A is greater than or equal to (h2 - c2), One may, through the use of 
Equation (9), select any number of ——- pairs (X,Y) which satisfy the 
ellipse determined by a chosen value of A(A2 h2 - c2). The value of A and 
any desired number of corresponding corrdinate pairs may readily be chosen 
in such a fashion that the quantity Ap/p(U2 + V2)/2 will be evaluated at a 
fairly representative number of points in the region of interest. All such 
points will lie on the ellipse defined by the selected value of A. If this pro- 
cess is repeated for several other choices of A it is possible to obtain a rea- 
sonably complete picture of the magnitudes of the pressure changes associ- 
ated with the pressure distribution. As will be seen by an inspection of 
Figure (2), this procedure has been followed in selecting nearly all of the 
coordinate points at which Ap/P(U2 + V2)/2 has been evaluated. 


— 


SUMMARY 


In the foregoing an attempt has been made to set forth a procedure which 
may be employed to obtain a rough approximation of the sea bottom pressure 
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distributions produced in the neighborhood of the upstream hull line by a sur- 
face vessel yawed more than 20°. The vessel is considered to move at low 
speed in a calm sea of moderate depth. To this end a treatment of the flow 
pattern around an ellipsoid model in an infinite nonviscous fluid half-space is 
given. A brief derivation of the most important formulae for the velocity and 
pressure distributions produced by the ellipsoid on the plane boundary is pre- 
sented. The mathematical complexities involved if one attempts to determine 
an analytical expression for the curves of equipressure change are pointed 
out and an alternative procedure for obtaining these curves is described. An 
application of the theory has been made to a special case. The results of cer- 
tain calculations made in this regard have been presented to illustrate the 
principal features of the pressure distributions which may be expected. 
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THE BOUNDARY LAYER DEVELOPMENT IN OPEN CHANNELS?# 


Discussion by Yoshiaki Iwasa, Alfred C. Ingersoll and Masashi Hom-ma 
Closure by J. W. Delleur 


YOSHIAKI IWASA.*—Prof. Delleur’s study of the growth of the turbulent 
boundary layer in open-channel flows as well as that by Dr. Halbronn is most 
outstanding. The formulation of the dynamic principles of boundary layer 
and of the behavior of open channel flows will be shown in a different way. 
Use will be made of the one-dimensional approach of classical hydraulics 
in the light of modern hydro-dynamics. 

The writer has studied for several years the problem of boundary layer 
growth associated with open channel flow; the results for a smooth channel 
are to be published in the Memoirs of the Fac. Eng., Kyoto University. This 
writers’ study primarily with the hydraulic behavior of boundary layer growth 
for initially rapid flows, whereas Prof. Delleur was concerned with the flow 
on a horizontal bed and with the hydraulic behavior of the main flow near 
critical regime. 

The writer wishes to show both theoretical and experimental results ob- 
tained from a study of boundary layer growth in open channel flow. 

This earlier theoretical study has been presented to the Annual Convention 
of JSCE in Kwansai District held at Kobe University, Kobe, Japan in 1955; the 
experimental research has been conducted at the Hydraulics Laboratory, 
Kyoto University, with the use of a very smooth lucite flume 36 ft. long and 
of slope varying from 0° to 30°. The main part of the experimentation was 
conducted by Messrs. H. Matsunami and S. Kinukawa in 1956 when the writer 
was at the Hydrodynamics Laboratory, Massachusetts Institute of Technology. 

Although our theoretical procedure is the same as that of Halbronn and 
Delleur, basic relationships of both main and boundary layer flows used by 
Prof. Delleur are hydrodynamically proved in what follows. 

Dividing the flow into two parts: main and boundary layer, the principle of 
conservation of momentum or the integration of the Eulerian equation of mo- 


tion for two dimensional main flow yields, under the assumption of hydro- 
Static distribution of pressure: 


where uy = velocity of main flow, x = distance along the channel bottom, h = 
fluid depth, g = acceleration of gravity, and @ = inclination angle of channel bed. 


a. Proc. Paper 1138, January, 1957, by J. W. Delleur. 


* Asst. Prof. of Hydraulics, Dept. of Civ. Eng., Kyoto University, Kyoto, 
Japan. 
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Integrating Eq. (1) with respect to x: to result in 


(2) 


hcos@ X sind 


2 
24 


This is the well known equation of Bernoulli and the basic equation of Prof. 
Delleur. 


The same principle used for boundary layer flow yields 


( 


= (sind - coset ) + U, d \udy (3) 


) 


where T, = boundary shear and 6 = boundary layer thickness. 

In order to derive the usual form of the boundary layer equation in terms 
of displacement and momentum thicknesses, Eq. (1) is multiplied by §, and 
then added to Eq. (3): 


d du, (4) 
P = fucy + U) dy - 


0 


With the aid of expressions of displacement and momentum thicknesses, Eq. 
(4) becomes 


G d9 du, 
dx * Wax (296) 


Thus, basic relations used by Prof. Delleur are shown to be valid with the one 
dimensional procedure of approach from hydrodynamics. 

The mutual relationship between the boundary layer approach and the usual 
one dimensional approach for mean flow, which was not explained by Prof. 
Delleur, can be established as follows: 

Integrating Eq. (1) with respect to y throughout the whole region of main 
flow, and adding this relation to Eq. (4) yields 


= ghsing - ghcos 4 ah (6) 
ax 


where 8 = momentum correction factor. Therefore, with the use of normal 
depth h, and critical depth h,, we can get 


dh h 
——= . (7) 
ax an @ A. 


This is the hydraulic equation derived usually by the one dimensional approach 
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for mean flow. It means that the introduction of the boundary layer concept to 
open channel flows does not require any new knowledge for the hydraulic 
equation in classical hydraulics. However, the transitional behavior of flow 
near an entrance to a channel is more precisely expressed by solving simul- 
taneously three basic equations for main flow, boundary layer, and constant 
discharge together. 

At the entrance to the channel there is a reach where the assumption of 
flow without curvature is not valid as this hypothesis makes the velocity near 
the bottom faster than that in the upper flow. Therefore, to make the actual 
flow consistent with the mathematical assumption at the initiation of growth 
used by Prof. Delleur and the writer, a lucite guide vane was set up for our 
experiments to get a uniform velocity distribution at the entrance zone. Al- 
though the influence of setting the guide vane on flow characteristics is ex- 
pected, the writer will only describe the following experimental results of 
velocity profile and its related skin friction law. 


1. For a smooth boundary, the shape of the velocity profile of the turbulent 
layer is better approximated by the power type of law than the logarithmic 
law as Dr. Bauer* and Prof. Delleur concluded. Although the behavior of H 
expressed by 6,/ indicates a quite different character from the usual con- 
clusion obtained by the unconfined theory, as seen in Fig. 1, it rapidly tends 
to a certain definite value determined by flow characteristics. 


2. The local skin friction coefficient, Ce = 2 To/R u2, has slightly lower 
values than those derived from the curves of Blasius and Dr. Bauer (Fig. 2). 


3. The Karman constant in the logarithmic velocity law decreases rapidly 
with the increase of Reynolds numbers and converges to a constant value of 
0.40, as seen in Fig. 3. 


Although the conclusion may be largely influenced by the guide vane, we 
can show the same behavior for steady uniform flow after the fluid has 
travelled downstream; this was shown by Dr. Keulegan and Prof. Powell. 

The mathematical behavior in terms of a power law derived by three basic 
equations indicates that the expression has a singular point very near the 
origin of growth. However, as an approximation, the development of boundary 
layer will be expressed in a form of 


_ Jr 


where A = ratio of skin friction coefficient to mth power of Reynolds numbers, 
v' = dimensionless parameter of kinematic viscosity and ¥ and o = the same 
parameters as Prof. Delleur used. For the case of the Blasius 7th law, 


29 


Fig. 4 describes the growth of boundary layer calculated by several 
approaches. It is seen that for a smooth boundary the blasius 7th law indi- 
cates a close approximation of actual development. 
Because we are concerned with the practical application to channel design, 
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the Manning roughness and the location of critical point where the boundary 
layer intersects the free surface should be considered. Fig. 5 indicates the 
behavior of the Manning roughness with an increase of distance and it is un- 
derstood that the Manning roughness has a constant value determined by 
channel characteristics, while it has a large value at the earlier stage of 
development. Therefore, we may conclude that the channel is adequately de - 
signed by the usual method of one dimensional approach, as a close approxi- 
mation. 

The determination of the location of the critical point is also of importance 
for air-entrained flow and laboratory investigation. Fig. 6 describes the lo- 
cation of the critical point of flows over a concrete shute under the assump- 
tion of Manning roughness n = 0.013 and E = (q/2)2/3. The plottings are ob- 
tained by field observation of Lovely and Michels. 

It is hoped that much of the true nature of fluid motion still unsolved will 
be soundly clarified through the physical knowledge obtained by the boundary 
layer theory. However, a clear formulation of the boundary layer process is 
quite difficult because of the great difficulty with the mathematics and with 
the systematic experimentation. 


ALFRED C. INGERSOLL,* M. ASCE.—One of the new horizons opened by 
this paper is the application of Mr. Delleur’s boundary layer theory to the 
determination of bed shear stresses in open channels with movable beds. 
Knowing only the channel geometry, discharge, surface velocity and water 
temperature, it is possible to compute values of Tp without having any mea- 
surement of the slope of the energy gradient. This attack is especially help- 
ful in the study of bed shears in settling tanks where the velocities are so low 
that no determination of slope can be made. 

The writer has had occasion to test out this application in a glass-walled 
flume 1.27 ft. wide by 14 ft. long, at depths varying from 0.610 to 0.623 ft., 
the floor being partially covered with particles of gilsonite ranging from 
0.0232 to 0.0328 in. in diameter. Surface velocities, uj, were measured by 
timing confetti streamers. For the purpose of this computation Eq. (17) is 
written in the more convenient form, 


2, (1) 


* 


uy 


in which V is the mean velocity and R is the hydraulic radius. The bed shear 
stress, Ty, is computed from Eq. (7) with m = 1/8. The measurements re- 
ported below were all taken at a point 8 ft. downstream from the inlet to the 
channel, a perforated plate baffle with about 10 per cent free area. For com- 
parison purposes, values of T, are also computed from the relation, To = 
{V2/8, where f is determined from the Moody chart of friction factor for 
pipes, (Ref. A) with Reynolds number equal to 4RV/y, and relative roughness 
given by the size of gilsonite particles relative to 4R, or 0.0018. 


* Associate Prof. of Civ. Eng., Calif. Inst. of Technology, Pasadena, Calif. 
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Mean Surface Displacement 3ed Shear Bed Shear by 
Velocity Velocity Thickness by Eq.(7) Friction Factor 
* 
fps fps g. psf psf 


i12 0.122 0.0270 0.925x10°* 0.924x1074 
0.116 6. 127 0.0272 0.985 0.991 
0.122 0.133 0.0258 1.093 1.094 
0.128 0.140 0.0268 1.180 1.199 


The remarkable agreement between the values of bed shear computed by 
these two methods is perhaps not surprising inasmuch as each method stems 
from the application of pipe flow formulas to open channels. The boundary 
layer approach, however, utilizes the Blasius formula and the concept of the 
virtual channel, while the friction factor method is based on the Colebrook- 
White function for flow in rough pipes. 

The concept of the virtual channel is interesting and helpful. If the equa- 
tion for the shaded corner areas if y'z' = 9.3126 ,2 as indicated in Fig. 1, 
however, the writer believes that the last term of Eq. (12) would come to 
-23 6,“ instead of the 0.56 ,2 shown. Even with this value, however, the term 
would still be less than 5 per cent of A, which might justify its omission for 
many types of computations. 

With reference to curve II in Fig. 3, the writer is confused by the state- 
ment that “rj is the value of r which satisfies Eq. (31):” Inasmuch as Eq. (31) 
is simply the equation of curve II, then all values of r will satisfy the equation 
for corresponding values of g. The values of rg, rj, rg, and rg in Fig. 3 all 
seem to be related to a value of o = 0.75. Whether this is simply a repre- 
sentative value of o or whether it has some special significance is not imme- 
diately clear. 

In Table I it is not clear what the difference is between the first and third 
rows. In both cases, presumably, the values of 6, have been computed from 
velocity profiles, but the values in the first row seem to be just 1.74 times 
those of the third row in each instance. 

The writer finds the concept of the “vertical tangent” referred to in the 
section on general properties of the boundary layer development to be a very 
difficult one. Regardless of what the theory may predict, the physical concept 
of a boundary layer above a horizontal bed having a vertical tangent seems 
indeed remote from reality. 

The author makes considerable use of H as a length parameter. As it is 
tied to the specific energy, Eg, it must vary along the length of a horizontal 
channel unless the reference value of Eg is fixed at one particular point. 
Similarly the author refers to “uniform flow” as that occurring in the horizon- 
tal channel after the boundary layer has pierced the surface. As it is com- 
monly understood, uniform flow cannot occur at any time in a horizontal chan- 
nel because the frictional losses cannot be equated to the drop in bed level. 
Clarification on these points will be helpful. 

In the analysis of experiments the author remarks that ¢ is not a free 
parameter when Eg and L are given. Inspection of the summary data reveals 
that Series B and C have very closely the same energy (0.259 against 0.2505) 
while their values of o vary widely (from 0.50 to 0.855). Considering Series 
A and C, on the other hand, reveals that the values of 0 are almost the same 
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while the specific energies differ considerably. Is it not true that the tank 
length, L, was the same for all of these tests? If it was, as would appear to 
have been the case, how may these anomolous relations be accounted for? 

This paper is a signal contribution to the understanding of open channel 
fluid mechanics, yielding as it does a means for computing the boundary layer 
development in rectangular channels. The writer hopes that the studies may 
be continued and extended to channels of other shapes, and especially to the 
very tranquil flows encountered in settling tanks. While still in the turbulent 
range, such flows would be characterized by values of the discharge parame- 
ter o on the order of 0.01, far below the range of values presented here. 

The writer wishes to acknowledge with thanks the assistance of Mr. A. R. 


Porush and Mr. H. J. Hansen in making the bed shear tests reported herein- 
before. 
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MASASHI HOM-MA.!—The author’s treatment of the problem is very in- 
teresting and valuable for the promotion of further research. The writer 


wishes to consider these phenomena and to suggest research on related 
problems. 


The author’s fundamental equations for the main flow may be written as, 


23 s 


=4,1h-d,) 


If we limit the problem to two-dimensional motion, these are the equations 
for the flow of an ideal fluid over a solid bottom, the height of which coincides 
with the displacement thickness in each section. Therefore, when we consider 
a flow which changes from subcritical to supercritical on this bottom, the 
control section must be the section where the value of Es- 6, is minimum. 
This fact shows that in such flow the displacement thickness must take its 
maximum at the control section. 

Actually, the discharge of ideal fluid through the depth h- 6, is not equal to 
¢ but is equal to he 4 (see figure 2). In the boundary layer there must occur 
the energy consumption which will cause the inflow of momentum from the 
main flow. (This must be the inflow of a certain mass of fluid from the main 
flow into the boundary layer). Thus the value of A q, will increase with the 
distance x from the inlet section. 

In the subcritical flow, the decreasing of 6, with the distance x is accom- 
panied with the increasing of the depth h-6, and the decreasing of uy (see 
Fig. 1). This means that the increasing of 4 g is impossible and also the de- 
creasing of 6, can’t occur in the subcritical flow. 

On the contrary, it is understood from Fig. 1 that the decreasing of § , can 
be accompanied with the increasing of Ag in the supercritical flow. Now the 


writer finds that the author’s idea is reasonable and it will lead us to valuable 
results. 


1. Prof. Dept. of Civ. Eng., Univ. of Tokyo, Tokyo, Japan. 
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Figure 2 
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But a doubt still remains. If we consider the case when the length ({) of 
the dam crest is very long, the author’s result shows that 6, continues to de- 
crease until the amount of uj becomes extremely large. Such velocity seems 
not to take place. Also, in an inclined open channel, an uniform flow state of 
ideal fluid is impossible, because the ideal fluid must be accelerated constant - 
ly in the channel. Rarely can we observe a fairly long extent of undisturbed 
surface in a steep-sloped flow, but when the flow becomes uniform, the 
boundary layer must extend itself to the free surface. It seems to the writer 
that after a long run of flow the boundary layer grows with increasing 47 
(discharge through the depth 6 ,) even in supercritical flow. Here we will 
look back the theoretical analysis of the problem. In the writer’s opinion, for 
the case of a long channel, a correction must be added to the pressure term 
in the equation (4) of the original paper, because the water depth above the 
boundary layer is not constant in this problem. 


J. W. DELLEUR,* A.M. ASCE.—The discussions by Messrs. Blench, 
Ingersoll, Iwasa and Hom-ma were most stimulating. Professor Blench 
makes a number of inquiries from the point of view of the practicing engineer, 
whereas, Messrs. Ingersoll and Iwasa present some additional experimental 
data, and, as Professor Hom-ma, they discuss the problem from a fluid me- 
chanics point of view. 

Consider first Professor Blench’s questions. A fully developed flow does 
not occur in the throat of a broad crested weir under usual conditions. The 
flow is critical, and the boundary layer reaches a maximum at a point near 
the brink. A fully developed flow can occur, however, in a horizontal open 
channel if it is long enough. The fully developed flow will start where the 
boundary layer reaches the surface. At this point, the value of the independent 
variable r is given by eq. (33), which replaced in (27) with the proper value of 
the flow parameter g, will yield the corresponding value of m. The value of 
m= 0.13 has no particular significance. 

Regarding Professor Blench’s second and third questions, eq. (15) defines 
the energy Eg of the potential flow, and H; it does not compel critical flow. 
The theory developed does not require that the flow be very nearly critical. 

It is solely pointed out that the boundary layer development will reach a maxi- 
mum when the flow is close to critical, as given by (34). 

In his questions (iv) and (v) Professor Blench challenges the use of the 
boundary layer momentum eq. (4) for turbulent flows, because of the neglect 
of the turbulent shear-stress. Although the justifications for so doing were 
already briefly mentioned in the paper immediately after eq. (4), it may be 
added here that Ross(5) and Robertson(6) have integrated the complete Rey- 
nolds equations, and have obtained four corrective terms due to anisotropy of 
the turbulence normal stresses, to the second derivative of the shear stress, 
to the deviation of the mean flow from Prandtl’s assumptions, and to the vis- 
cosity. Essentially these corrections are all connected with the fact that the 
pressure is not exactly constant across the boundary layer. With the first 
correction, which is the predominant one, the momentum equation reads 


J 
5; (4 ~ ) dy 


* Associate Prof. of Hydr. Eng., School of Civ. Eng., Purdue University, 
Lafayette, Ind. 
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However, past experiments have shown that this correction may be neglected 
for zero or small pressure gradients. 

Finally, it is not agreed, as Professor Blench states in his sixth question, 
that 0 is H3/2 divided by itself. o has the form of a Froude number, and 
perhaps, isan unconventional one. However, there are many accepted forms 
of the Reynolds number, in which the length term may be a pipe diameter, a 
displacement thickness, a grain size, or the velocity term may be /To/?p, 
etc. Similarly, there is no reason why different forms of the Froude number 
cannot be constructed, as long as its form and non-dimentionality are re- 
spected. 

The discussion by Dr. Ingersoll is most interesting, in particular because 
of his novel use of the boundary layer approach to compute the bed shear 
stress in open channels with movable beds. He computes the boundary layer 
displacement thickness from 


which is the same as (26) in which the hydraulic radius, R, replaces the depth. 
His comparison of the bed shear stress obtained by Blasius formula, eq. (7), 
and the Colebrook-White function for pipe flow is noteworthy. 

One of the difficulties encountered in the computation of the virtual channel 
area is that of evaluating the zone of mutual interference of the displacement 
thicknesses in the corners of the channel. In view of the question asked by 
Dr. Ingersoll, the detail of the computation which was omitted in the paper 
for the sake of brevity, is given in what follows. Loitsianskii and Bolshakov(9) 
modified the 1/7th power law for the velocity distribution for the zone of in- 
terference of two boundary layers near the intersection of two plates, (fig. 
D-1), to 


[ (D-1) 
a(x, 8) 
and Blasius formula to 
tT 0.0225 Vy (D-2) 
eu? sine 


where x is measured from the leading edge of the plates forming the corner 
in the direction of the flow, and @ is the angle of the corner. They show that 
the equation of the edge of the boundary layer is 


a(x,@) (D-3) 
where 
a(x,6) K(9) (x) (D-4) 
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in which K(@) is a certain non-dimensional constant parametrically dependent 
on the angle 6 , and 4 o is the boundary layer thickness in the undistrubed re- 
gion, away from the corner. The constant K(@) has been determined by 
Loitsianskii and Bolshakov,(9) and for @ = 90°, K(@) = 1.164. Now, letting 
Zo(x) be the outside limit of the zone of interference, measured along the side 
or the bottom of the channel from the corner, at station x, as shown in fig. 
(D-1). For y = Zo(x), Z' = 6 (x), equations (D-3) and (D-4) yield 


le (x) a (x,9) K (6) 5 *(x) 
§ (x) y 


Since K (90°) = 1.164 and since for the seventh power law 6 = 86 ,, 
(x) = 8x 1164 (x) = 9.312 (x) 


Referring to figure (D-1) and calling A7, the total area OAB'C' DO, Ay, the 
“lost” area OABCDO, and Ay the net area BB'CC'B, it follows from the 
virtual channel concept that Aw = Ay - Ay, where Aw represents the area 
through which an ideal fluid flows with the same discharge as the real fluid 
through the total area. Thus 


why (D-6) 
uw 


where Q" is the discharge through Ay. The curved boundary B'C"' of the 
zone of interference is defined by (D-3) which becomes 


= a(x,6) = K(@) 8*(x) = 274-496 


Computing the area Ar: 


dz'd 85. B08 


72. 
A, 9.312 + | 


° 


where 


a; 74.496 


and using the modified 7th power law (D-1) in (D-6) 


74-496 


dztdy’} 


74.490 


ge: 
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[z+ In (1.164)] = 18.874 


The cross section of the imaginary channel is thus 


Az bh - (b- 18.6243) S¢ -2 (h-9.312 18.874 


or 


A= bh «be ~2hd« -0,500 3% (D-7) 


With reference to the curves of fig. 3, the values of rp, rj, rg, and rg are 
the values of r which satisfy equations (30), (31) and (35) for a given value of 
0. For each discharge there is a particular value of o for which there is a 
value of ro, ry, r2 andrg. These values are convenient for numerical com- 
putations, because they indicate the range of values that the independent vari- 
able r may take. The most reliable part of the curve of the growth of the 
boundary layer is that between its origin and its maximum-—if it exits—that 
is forrg< r<rj,. It was stated after eq. (35), and is emphasized again that 
for r > ry, the curve of growth of the boundary layer loses its accuracy, be- 
cause of the curvature of the flow, which becomes important in that reach, 
was not considered in the derivation. As r approaches rg, the theoretical 
curve ceases to represent accurately the physical phenomena. The plot of 
equations (27) and (28) see fig. 2) shows a point of vertical tangent, for r = ro, 
which, as Dr. Ingersoll points out, is “remote from reality.” Nevertheless, 
the important conclusion, which is confirmed by the experiments, is that be- 
yond the maximum, for which r = rj, the boundary layer decreases. 

The length parameter H is tied to the specific energy Eg, (see eq. 15). For 
the purpose of the computations Eg was measured at the beginning of the tur- 
bulent boundary layer. In the experiments the value of Es for series B and C 
are very close, but the value of the discharge parameter o varies consider- 
ably. In Series B the flow was controlled by an adjustable weir at the down- 
stream end of the channels, whereas in Series A the channel ended in a free 
fall. In all cases the length of the test flume was the same. The effect of the 
length of channel will be analyzed to show how Eg, ¢@, and L are related, and 
to show that there is no anomalous relation, as asked by Dr. Ingersoll. 

It was stated that the equations (19) and (24) governing the development of 
the boundary layer do not contain explicitly the length L of the channel, be- 
cause the energy Eg, the discharge parameter o and the length of the channel 
are not independent. If Eg and L are given, o may be computed. Two cases 
will be discussed, depending upon whether or not the maximum of the boundary 
layer exists (see eq. 34). Consider first the case of a broad crested weir, 
where the maximum exists. The length d between the origin and the point of 
maximum development may be obtained by integrating (28) between rg and rj, 


assuming c = constant, and making use of (30) and (31) for the lower and upper 
limits respectively: 
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The dotted line in fig. 2 shows the locus of the maxima computed by the above 
equation, referred to as the approximate formula because of the approxima- 
tion involved in taking c = constant. It should be noted here that c is very un- 
sensitive to changes in r and thus the approximation is a good one. On the 
other hand, d may be computed independently, by evaluating L - d, or the dis- 
tance between the point of maximum boundary layer or point of critical flow 
and the free overfall. For this purpose the criteria developed by Rouse(A) 
may be used. If Eg and L are given, H may be computed from (15), L - dis 
estimated from the criteria, d is thus known; (D-8) may be solved for ro, 
which, introduced in (30) will yield the value of 0. This reasoning may be 
used in developing a simple way of computing the discharge coefficient for 
flow over a broad crested weir.(16) 

A similar argument may be developed when the boundary layer pierces the 
free surface if the length from the origin of the layer to the point of fully 
developed flow and the energy Eg are known, it may be shown that the dis- 
charge parameter o may be evaluated. 

Professor Iwasa’s discussion is commended for rederiving the boundary 
layer momentum equaticn (4) from the one-dimensional analysis of open chan- 
nel flow, and for tying the boundary layer principles to the classical Bresse 
formula (7 of Iwasa’s discussion). Attention should be called to the unfortu- 
nate change of notation in Iwasa’s discussion, in particular in eq. (8) of the 
discussion, which is clarified in the table below: 


Delleur’s Iwasa’s 
Item Notation Notation 
shape arameter H 
critical depth H he 
6 * 
m no symbol 
Power of Lbs in (7) no symbol m 
- 1 in eq. (7) 
4 
width 
critical depth no symbol 
skin friction coefficient no symbol 7 


power of Reynolds No. in (7) 


Among the experimental data presented by Mr. Iwasa, fig. 4 of Iwasa’s 
discussion is of special interest because it compares the values of the de- 
velopment obtained by experiment and analytically by means of Blasius 7th 
power law. 

It is not agreed, as stated by Prof. Hom-ma, that the discharge through 
the depth h- §, is 47 44; but simply q as derived in (10). It is true, 
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however, that there is a flow from the main flow into the boundary layer, and 
this is considered in the derivation(3) of the boundary layer momentum 
equation (4). 


REFERENCE 


A. Rouse, H., “Verteilung der hydraulischen Energie bei einem lotrechten 
Abstrus, Oldenburg, 1933. 
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Fig. D-1. Corner of Channel 
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A PRESSURE LINE CONCEPT FOR INELASTIC BENDING* 
Closure by Frank Baron 


FRANK BARON, ! M. ASCE.—Thanks are given to Mr. Hribar for his dis- 
cussion and for his pertinent questions in regard to the pressure line concept 
for inelastic bending. 

Mr. Hribar agrees that the pictorial procedure can be useful in quickly 
estimating the effects of plasticity on the behavior of structural elements sub- 
jected to axial and flexural loads. He also indicates that an estimate is fre- 
quently all that is desired in studies of the inelastic behavior of structural 
frameworks. The writer wishes to add that the degree of accuracy desired 
can depend on the needs and demands of an analyst. It is observed that the 
needs of an analyst are not always the same and that they do not always war- 
rant “exact” solutions to problems. This is particularly true in the various 
stages of a design. 

The statement is made by Mr. Hribar that the procedure results in only an 
estimate unless an analyst is experienced in the method and has a keen sense 
of judgment in balancing statics pictorially. Fortunately, little experience and 
almost no judgment are required in satisfying the requirements of statics. 
The only judgment required is to observe that the correction diagram as de- 
fined in the paper vanishes or balances. The diagram is balanced when the 
positive areas equal the negative areas and the centroid of the positive areas 
coincides with the centroid of the negative areas. This requirement is so 
simple that in most cases it can be judged by eye. However, as an aid to 
judgment, a numerical check can be made of any given estimate. 

Several important attributes of the procedure need to be emphasized. In 
the theory of plasticity, the requirements of statics, geometry, and properties 
of materials are immediately apparent when stated pictorially and when com- 
pared to the results of the elementary theory of mechanics. In this way, 
every relationship can be observed including the relationship of the theory of 
plasticity to the elementary theory of mechanics. 

Concerning the specific questions raised by Mr. Hribar, not all members 
are made of steel with idealized stress-strain curves nor are all cross- 
sections of members sufficiently simple to lend themselves to algebraic 
formulation. However, if an algebraic formulation is desired for a given 
stress-strain curve and a given shape of cross-section, the concept presented 
in the paper can be used with the same ease as present methods. At times, it 
can even be used with greater ease. For example, consider a parabolic 
stress-strain curve and a simple rectangular cross-section. Now determine 
by means of the pressure line concept the value of the reduced modulus, Ep, 


a. Proc. Paper 1157, January, 1957, by Frank Baron. 
1. Prof. of Civ. Eng., Univ. of California, Berkeley, Calif. 
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for any combination of axial load and moment about a principal axis of the 
section. The answer is obvious when the properties of a parabola are 
recognized. 

As to the determination of ultimate loads, the same statements apply here 
as above. Mr. Hribar agrees that as the cross-section or the stress-strain 
curve becomes more unwieldy the pressure line concept becomes more 
powerful. No further comments are required by the author. 

In conclusion, the writer is grateful to Mr. Hribar for the opportunity to 
emphasize certain aspects of the pressure line concept. 
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DEFLECTIONS OF STRUCTURES IN THE INELASTIC RANGE4 
Discussion by John H. Percy 


JOHN H. PERCY.!—Mr. Gerstle’s paper may be considered in two parts, 
namely the computation of suitable coefficients to describe the behavior of 
elements of a structure and the use of charts showing these coefficients to 
determine deflections in any particular case. The second part holds no mat- 
ter what shape the curves take (provided always that the coefficients are 
Single valued). The treatment would therefore be more general if the charts 
could readily be drawn for any section shape and any stress-strain relation- 
ship. 

If @ is the curvature of a member subjected to a bending moment M the 
relationship between the two will, in general, have the form shown nondimen- 
sionally in figure (a). 


For figure 1 of the paper, taking the origin at A and measuring X along the 
beam, geometry gives 


L 
| pade Oa = | pdx. 


From these it follows, using 
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a. Proc. Paper 1290, January, 1957, by Kurt H. Gerstle. 
1. Lecturer in Civ. Eng., Univ. of Auckland, New Zealand. 
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where Cy, Co are the coefficients shown in figures 2 and 3 of the paper. 

Now, if gMyp is any bending moment and if, in figure (a), Pg is the area 
and Qg the first moment of area about the @ axis of the region OAB it follows 
immediately that 


h(Q,-Q,) —( FP, 
(h—-k)* 
Qn- Qe 
C. 


A similar analysis for figure 5 of the paper does not lead to such straight- 
forward results for 9 AC and 6 AC but it can be shown that 


a M n 
My, _'_lerg Alm.) [erp 
Sac were C,= Moe 2(n-k)} Myp 
2 n 
C, EI where C,= 2th-ky} 
Myp 
M,L Ale! 
=C, t= $32 
war C= 2fn-k) Mor 


(It is worth noting the rather surprising result that 6c, for this case and 
6AB for figure 1 differ only by a constant factor 5 for any bending moment- 
curvature relationship. The curves in figures 3 and 5 of the paper are there- 
fore essentially the same.) 

The coefficients of the paper may thus be calculated without difficulty for 
any bending moment-curvature relationship. 

This is best found by direct experiment for any particular material and 
section shape, since the determination of the stress-strain relationship with 
sufficient accuracy for use in calculating flexural deformations is a precision 
laboratory task requiring considerable care (c.f. Roderick and Phillipps). (i) 
If necessary however the bending moment-curvature relationship for any sec- 
tion may be derived from the stress-strain relationship by methods similar to 
those used by Hrennikoff, (ii) Phillips, (iii) Roderick, (¥) and others. 
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Fig. a. Relationship between Bending Moment and Curvature. 
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THE LATERAL RIGIDITY OF SUSPENSION BRIDGES# 


Discussion by Arne Selberg 


ARNE SELBERG.*—The paper demonstrates that the uniform distribution 
method is of little value for statical calculations. A short investigation will 
demonstrate that, for example, a triangular distribution will give a better re- 
sult. However, there is little reason for using any such approximate 
solutions. 

Theimer,(1) in 1941, has given the same solution for statical computations 
as given by Silverman in his Elastic Distribution Method. Theimer also used 
dimensionless parameters in his analysis. However, the solution as given by 
Silverman will be somewhat more convenient in use. 

The calculated example demonstrate a slow convergence for the series 
when calculating moments. 

Introducing (2,3) 


will improve the convergence. 

The reason is that the first series is integrated twice when finding the mo- 
ment, the second series has to be derivated twice. 

In his investigations on the oscillations in the horizontal plane. Silverman 
gives what appears to be an erroneous equation (7), and all results obtained 
from this equation are erroneous as well. It is easily demonstrated on a 
model that his results are incorrect. 

With the notations of Silverman, the equations for any oscillations will be: 


= -Z - 


a. Proc. Paper 1292, July, 1957, by I. K. Silverman. 


* Prof. of Theory of Structures, Norwegian Technical Univ., Trondheim, 
Norway. 
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Mg, and mc being mass of suspended structure and cables respectively 


Table III demonstrate clearly that mg and m¢ both are of importance 
From these equations it is impossible to get eq. (7) given by Silverman 


Introducing 


pli = SA, sin wet 


and 
U-Sinwe sin wel 
we get 
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These equations yield 


The coefficients A, are given by the equation 


2,? £ = 2 ty $14 


o 


The easiest solution is obtained by requiring that this equation shall be 
satisfied for distinct values of x. A somewhat better solution is given by re- 


quiring: 


Both methods give a set of linear equations for Ap, and w? is given by |D! 


= 0. 
With only one term 
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we have with the two methods: 


Method 1. 
quiring (Y% = at the middle of the span. 
Method 2. 
? 


The difference in results of the two methods are insignificant for most 
bridges. 

Both equations usually give two values of wy e.g. cable and truss oscillat- 
ing in same or opposite direction, only the first one is of special interest to bs: ‘ 
bridge designers. 

As the horizontal oscillations of a suspension bridge is of secondary im- q 
portance for the designer it will be little reason for using more than one or = 
two terms in an actual calculation. a 

For Tacoma Narrows (1940) the 2 formulas give 


fog 
w= 53) and W=4%2 (20.6) 


or 


"4420 /) fomin (204) 


respectively. 
A solution with two terms Aj; Ag give: 


N, = 41 min™! (22.9) and N, = 16.1 min”! (36.5) 


The difference in Ny is insignificant. The number in parentheses is for : 
cable and truss oscillating in opposite direction. Silverman give Ny = 9,9 a 
min~! for this bridge. 
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ANALYSIS OF CONTINUOUS BEAMS4 


Discussion by S. J. Medwadowski 


S. J. MEDWADOWSKI, ! A.M. ASCE.—Two principal advantages, achieved 
at a presumably slight loss of accuracy, are claimed by the author for the 
use of the Fourier series technique in the analysis of problems of continuous 
beams of constant bending rigidity. These are: 


a. the possibility of representing deflection, and related quantities, by means 
of a single expression valid for 0< x < ¢; 


b. the numerical work being reduced to a minimum. 


The first of these advantages results also by employing classical methods 
leading to a closed form solution if use is made of the Heaviside unit func- 
tion.2 The closed form solution yields exact (within the assumptions of the 
theory) results; in addition it possesses the advantage of yielding directly 
closed form expressions for bending moments and shears, as well as slope, 
by formally differentiating the expression for deflection. Because of poor 
convergence of the series involved these expressions have to be found by a 
separate, indirect procedure if the Fourier series representation is used, as 
pointed out by the author on p. 7 of the paper. 

Any saving in the numerical work involved in calculating the unknown re- 
actions Rj must occur in the calculation of the coefficients of Eqs. (15), in 
particular in calculating the coefficients Djj, defined in Eqs. (17). Written in 
a somewhat more compact notation the definition is: 


It is easy to show that the closed form solution leads to the following, 
essentially equivalent, definition: 


= ms (1 (2 = (i j (B) 


Ww 


a. Proc. Paper 1399, October, 1957, by Seng-Lip Lee. 

1. Head, Structural Dept., John A. Carollo, Cons. Engrs., Berkeley, Calif. 
2. A popular yet clear and concise presentation is given in “The Application 
of Heaviside’s Step-Function to Beam Problems,” by John E. Goldberg, 
Proc. ASCE, Separate No. 202, 1953. 
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In both cases Djj = Djj. It is readily seen that the numerical work of calcu- 
lating Djj from Eq. (A) involves eleven operations (if only two terms of the 

series are retained), while the use of Eq. (B) results in six operations only 
and yields an exact value of the coefficient. 

Thus it would appear that the advantages claimed by the author for the 
Fourier series representation are in fact illusory. 

The system of Eqs. (15) of the paper is usually derived directly from Eq. 
(11) by substituting in it the prescribed deformations y(ajf) = - dj. It would 
be of interest to learn the author’s reasons for making use of Castigliano’s 
theorem instead, a somewhat lengthier procedure. 

It is well known that a formal application of the Fourier series expansion 
technique often yields solution of a physical problem, even though some of 
the series involved in the process may be divergent. Bending of continuous 
beams is such a problem. However, the author’s statement of p. 3 of the 
paper, to the effect that a function may be “approximated” by a Fourier 
series, is incorrect. 

The application of the Fourier series representation to problems of static 
bending of continuous beams may be considered a particular case of the 
problem of dynamic bending, treated by Saibel and D’Appolonia.3 It is hoped 
that a more complete list of references will be provided by the author in his 
closing discussion. 


3. “Forced Vibrations of Continuous Beams,” by Edward Saibel and Elio 
D’Appolonia, Proc. ASCE, Separate No. 100, 1951. 
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